text:
# ## [Aeroelastic Analysis of a Typical Section](@id section-stability)
#
# In this example, we perform a two-dimensional aeroelastic analysis for a two degree of
# freedom typical section model.
#
# ![](../assets/section-drawing.svg)
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`section-stability.ipynb`](@__NBVIEWER_ROOT_URL__/examples/section-stability.ipynb).
#-
#
# The equations of motion for this model are
# ```math
# m \left(\ddot{h}+b x_\theta \ddot{\theta} \right) + k_h h = -L \\
# I_\theta \ddot{\theta} + m b x_\theta \ddot{h} + k_\theta = M
# ```
# where ``b`` is the semichord length, ``k_h`` is the linear spring constant, ``k_\theta``
# is the torsional spring constant, ``m`` is the mass per unit span, ``x_\theta`` is the
# distance to the center of mass from the reference point, ``I_θ`` is the moment of
# inertia about the reference point, ``L`` is the lift per unit span, and ``M`` is the
# moment per unit span about the reference point.  To learn more about how this model is
# implemented in this package, see [`Section`](@ref).
#
# This package has a number of pre-implemented aerodynamic models which may be used to
# model the aerodynamics of the typical section model.  These models include a steady-state
# thin airfoil theory model (see [`Steady`](@ref)), a quasi-steady thin airfoil theory
# model (see [`QuasiSteady`](@ref)), an unsteady aerodynamic model based on Wagner's
# function (see [`Wagner`](@ref)), and Peters' finite state aerodynamic model
# (see [`Peters`](@ref)).  We will perform aeroelastic analyses using each of these models
# and compare the results.
#
# The non-dimensional parameters we will use are
# ```math
# a = -1/5 \quad e = -1/10  \\
# r^2 = \frac{I_P}{m b^2} \quad \sigma = \frac{\omega_h}{\omega_\theta} \\
# \mu = \frac{m}{\rho_\infty \pi b^2} \quad V = \frac{U}{b \omega_\theta}
# ```
# where ``a`` is the normalized distance from the semichord to the reference point, ``e``
# is the normalized distance from the semichord to the center of mass, and ``\omega_h`` and
# ``\omega_\theta`` are the uncoupled natural frequencies, defined as
# ```math
# \omega_h = \sqrt{\frac{k_h}{m}} \quad \omega_\theta = \sqrt{\frac{k_\theta}{I_P}}
# ```

using Aeroelasticity, DifferentialEquations, LinearAlgebra

## define non-dimensional parameters
V = range(1e-6, 3.1, length=1000) # = U/(b*ωθ) (reduced velocity)
a = -1/5 # reference point normalized location
e = -1/10 # center of mass normalized location
μ = 20 # = m/(ρ*pi*b^2) (mass ratio)
r2 = 6/25 # = Iθ/(m*b^2) (radius of gyration about P)
σ = 2/5 # = ωh/ωθ (natural frequency ratio)
xθ = e - a # distance from center of mass to reference point
a0 = 2*pi # lift curve slope
α0 = 0 # zero lift angle
cd0 = 0 # drag coefficient
cm0 = 0 # moment coefficient

## choose dimensional parameters
b = 1 # semichord
ρ = 1 # air density
ωθ = 1 # pitch natural frequency
c = 343 # air speed of sound

## calculate dimensionalized parameters
U = V*b*ωθ # freestrean velocity
m = μ*ρ*pi*b^2 # mass
Sθ = m*xθ*b # mass imbalance
Iθ = r2*m*b^2 # inertia
ωh = σ*ωθ # plunge natural frequency
kh = m*ωh^2 # plunge spring constant
kθ = Iθ*ωθ^2 # pitch spring constant

## define aerodynamic models
aerodynamic_models = (Steady(), QuasiSteady(), Wagner(), Peters{6}())

## initialize eigenvalue/eigenvector storage
λ = Vector{Matrix{ComplexF64}}(undef, length(aerodynamic_models))
Uλ = Vector{Array{ComplexF64,3}}(undef, length(aerodynamic_models))
Vλ = Vector{Array{ComplexF64,3}}(undef, length(aerodynamic_models))

## perform an analysis for each aerodynamic model
for imodel = 1:length(aerodynamic_models)

    ## define coupled model
    model = assemble_model(;
        aerodynamic_model = aerodynamic_models[imodel],
        structural_model = Section())

    ## define ODE function
    f = ODEFunction(model)

    ## eigenvalue/eigenvector storage
    nλ = number_of_states(model)
    λ[imodel] = zeros(ComplexF64, nλ, length(V))
    Uλ[imodel] = zeros(ComplexF64, nλ, nλ, length(V))
    Vλ[imodel] = zeros(ComplexF64, nλ, nλ, length(V))

    ## loop through each reduced frequency
    for i = 1:length(V)

        ## define aerodynamic parameters
        aerodynamic_parameters = (; a = a, b = b, a0 = a0, alpha0 = α0, cd0 = cd0, cm0 = cm0)

        ## define structural parameters
        structural_parameters = (; kh = kh, ktheta = kθ, m = m, Stheta = Sθ, Itheta = Iθ)

        ## define additional parameters
        additional_parameters = (; U = U[i], rho = ρ, c = c)

        ## define parameter vector
        p = assemble_parameters(model;
            aerodynamic_parameters = aerodynamic_parameters,
            structural_parameters = structural_parameters,
            additional_parameters = additional_parameters)

        ## define initial guess for equilibrium states
        x0 = assemble_states(model)

        ## find equilibrium point
        x = solve(SteadyStateProblem(f, x0, p))

        ## linearize about equilibrium point
        K, M = linearize(model, x, p)

        ## perform linear stability analysis
        λi, Uλi, Vλi = get_eigen(model, K, M)

        ## correlate eigenvalues
        if i > 1
            ## previous left eigenvector matrix
            Uλpi = Uλ[imodel][:,:,i-1]

            ## use correlation matrix to correlate eigenmodes
            perm, corruption = correlate_eigenmodes(Uλpi, M, Vλi)

            ## re-arrange eigenmodes
            λi = λi[perm]
            Uλi = Uλi[perm,:]
            Vλi = Vλi[:,perm]
        end

        ## save eigenvalues/eigenvectors
        λ[imodel][:,i] = λi
        Uλ[imodel][:,:,i] = Uλi
        Vλ[imodel][:,:,i] = Vλi
    end
end

#!jl nothing #hide

# We now plot the results for each aerodynamic model.

using Plots
pyplot()

sp1 = plot(
    title = "Non-Dimensional Frequency",
    xlim = (0,3.1),
    xtick = 0.0:0.5:3.0,
    xlabel = "\$ \\frac{U}{b \\omega_\\theta} \$",
    ylim = (0, 1.05),
    ytick = 0.0:0.2:1.0,
    ylabel = "\$ \\frac{\\Omega}{\\omega_\\theta} \$",
    titlefontsize = 10,
    guidefontsize = 10,
    legendfontsize = 8,
    tickfontsize = 9,
    foreground_color_legend = nothing,
    background_color_legend = nothing,
    minorgrid=false
    )

sp2 = plot(
    title = "Non-Dimensional Damping",
    xlim = (0,3.1),
    xtick = 0.0:0.5:3.0,
    xlabel = "\$ \\frac{U}{b \\omega_\\theta} \$",
    ylim = (-0.7, 0.605),
    ytick = -0.6:0.2:0.6,
    ylabel = "\$ \\frac{Γ}{\\omega_\\theta} \$",
    framestyle = :zerolines,
    titlefontsize = 10,
    guidefontsize = 10,
    legendfontsize = 8,
    tickfontsize = 9,
    legend = :topleft,
    foreground_color_legend = nothing,
    background_color_legend = nothing,
    minorgrid = false
    )

labels = ["Steady", "Quasi-Steady", "Wagner", "Peters (N=6)"]

for ia = 1:length(aerodynamic_models)

    plot!(sp1, V, imag.(λ[ia][1,:])/ωθ,
        label = labels[ia],
        color = ia,
        markersize = 1,
        markerstrokewidth = 0,
        )

    for i = 2:size(λ[ia], 1)
        plot!(sp1, V, imag.(λ[ia][i,:])/ωθ,
            label = "",
            color = ia,
            markersize = 1,
            markerstrokewidth = 0,
            )
    end

    plot!(sp2, V, real.(λ[ia][1,:])/ωθ,
        label = labels[ia],
        color = ia,
        markersize = 1,
        markerstrokewidth = 0,
        )

    for i = 2:size(λ[ia], 1)
        plot!(sp2, V, real.(λ[ia][i,:])/ωθ,
            label = "",
            color = ia,
            markersize = 1,
            markerstrokewidth = 0,
            )
    end
end

p1 = plot(sp1, sp2, layout = (1, 2), size = (800, 300))

#md savefig(p1, "../assets/section-stability.svg") #hide
#md nothing #hide

#md # ![](../assets/section-stability.svg)

# Using the `Wagner` or `Peters` aerodynamic models yields a flutter reduced velocity
# around 2.2, while the `Steady` and `QuasiSteady` aerodynamic models predict
# significantly lower flutter velocities.  The aerodynamic state variables of the
# `Wagner` and `Peters` models allows these models to capture the impact of vortex
# shedding on the lift and drag of the profile, therefore we can expect these models to
# yield more accurate results than the `Steady` and `QuasiSteady` models.
#
# The non-dimensional parameters we use for this example match those used by Hodges and
# Pierce in "Introduction to Structural Dynamics and Aeroelasticity".  Hodges and Pierce
# performed the analysis using a steady-state model and Peter's finite state model with
# six state variables.   The results presented here for the steady-state and Peters'
# finite state models match the results presented by Hodges and Pierce in "Introduction
# to Structural Dynamics and Aeroelasticity", which validates our implementation of
# these models.  Additionally, since the flutter speed predicted by the `Wagner` and
# `Peters` models match, we can be reasonably confident that the Wagner unsteady
# aerodynamic model is also implemented correctly.


text:
"""
    Fermi Quantum Chemistry Module

Authors: G. J. R. Aroeira, M. M. Davis, J. M. Turney, and H. F. Schaefer

GitHub: [Fermi.jl](https://github.com/FermiQC/Fermi.jl)
"""
module Fermi

using LinearAlgebra: Threads
"""
    Fermi.AbstractWavefunction

Abstract type common to all wave functions.

_struct tree:_

**AbstractWavefunction**  (Top level)
"""
abstract type AbstractWavefunction end

include("Backend/Options.jl")
include("Backend/Arrays.jl")
include("Backend/Contract.jl")
include("Backend/Output.jl")
include("Backend/Libcint.jl")
include("Backend/PhysicalConstants.jl")
include("Core/DIIS.jl")
include("Core/Molecule.jl")
include("Core/BasisSet.jl")
include("Core/Orbitals.jl")
include("Core/IntegralHelper.jl")
include("Methods/HartreeFock/HartreeFock.jl")
include("Methods/MollerPlesset/MollerPlesset.jl")
include("Methods/CoupledCluster/CoupledCluster.jl")

include("Tools/SinglePointEnergy.jl")
include("Tools/FiniteDifferences.jl")

end # module

ext: jl


text:
## -------------------------------------------------------------------
# Glc uptake
function plot_residual_Glc_study()

    p = plot(;
        xlabel = _textbf("experiment "),
        ylabel = _textbf("s_g", "/", "c_g"),
    )

    for Ed in [ChK.KayserData, ChN.NanchenData, ChF.FolsomData]

        src = nameof(Ed)
        marker = (9, source_markers[Ed])

        DAT = ChE.load_DAT(src)
        EXPS = DAT[:EXPS]

        @info("Doing", src)

        sGLC = Ed.val.(:sGLC, EXPS, 0.0)
        cGLC = Ed.val.(:cGLC, EXPS, 0.0)

        scatter!(p, EXPS, abs.(sGLC ./ cGLC);
            label = "",
            color = :black,
            marker,
        )
        # plot!(p, EXPS, abs.(sGLC ./ cGLC);
        #     label = "",
        #     color = :black,
        #     ls = :dash,
        #     alpha = 0.8,
        # )

    end

    # plot!(;legend = :bottomright)
    # sfig(ChE, p,
    #     "glc_residual_study", ".png"
    # )

    p
end


text:
### Quaternion type heavily based off of Base lib Complex.jl

"""
    Quaternion{T<:Real} <: Number
Quaternion number type with real and imaginary parts of type `T`.
`QuatF64` and `QuatF32` or `QuaternionF64` and `QuaternionF32` are aliases for
`Quaternion{Float64}` and `Quaternion{Float32}` respectively.
"""
struct Quaternion{T<:Real} <: Number
    re::T
    im::T
    jm::T
    km::T
end

### Constructors
Quaternion(a::Real, b::Real, c::Real, d::Real) = Quaternion(promote(a,b,c,d)...)
Quaternion(x::Real) = Quaternion(x, zero(x), zero(x), zero(x))
Quaternion(z::Complex) = Quaternion(real(z), imag(z), zero(real(z)), zero(real(z)))
Quaternion(q::Quaternion) = q

(::Type{Quaternion{T}})(x::Real) where {T} = Quaternion{T}(x, 0, 0, 0)
(::Type{Quaternion{T}})(z::Complex) where {T} = Quaternion{T}(z.re, z.im, 0, 0)
(::Type{Quaternion{T}})(q::Quaternion) where {T} = Quaternion{T}(q.re, q.im, q.jm, q.km)

### For representation as q = z + c*j for complex z and c
Quaternion(z1::Complex, z2::Complex) = Quaternion(real(z1), imag(z1), real(z2), imag(z2))

# Doesn't seem like a good design, fix this ?
Quaternion(v::Vector{<:Real}) = Quaternion(v[1], v[2], v[3], v[4])

"""
    jm
Quaternion imaginary unit.
"""
const jm = Quaternion(false,false,true,false)
"""
    km
Quaternion imaginary unit.
"""
const km = Quaternion(false,false,false,true)

const QuaternionF16 = Quaternion{Float16}
const QuaternionF32 = Quaternion{Float32}
const QuaternionF64 = Quaternion{Float64}
const QuatF16 = Quaternion{Float16}
const QuatF32 = Quaternion{Float32}
const QuatF64 = Quaternion{Float64}

convert(::Type{Quaternion{T}}, x::Real) where {T<:Real} = Quaternion{T}(x,0,0,0)
convert(::Type{Quaternion{T}}, z::Complex) where {T<:Real} = Quaternion{T}(real(z),imag(z),0,0)
convert(::Type{Quaternion{T}}, q::Quaternion) where {T<:Real} = Quaternion{T}(q.re, q.im, q.jm, q.km)

promote_rule(::Type{Quaternion{T}}, ::Type{S}) where {T<:Real,S<:Real} =
    Quaternion{promote_type(T,S)}
promote_rule(::Type{Quaternion{T}}, ::Type{Complex{S}}) where {T<:Real,S<:Real} =
    Quaternion{promote_type(T,S)}
promote_rule(::Type{Quaternion{T}}, ::Type{Quaternion{S}}) where {T<:Real,S<:Real} =
    Quaternion{promote_type(T,S)}

widen(::Type{Quaternion{T}}) where {T} = Quaternion{widen(T)}

bswap(q::Quaternion) = Quaternion(bswap(q.re), bswap(q.im), bswap(q.jm), bswap(q.km))

### Redefine/overload methods for Complex
real(q::Quaternion) = q.re
imag(q::Quaternion) = q.im

"""
    jmag(q)
Return the j-imaginary part of the quaternion number `q`.
# Example
```jldoctest
julia> jmag(1 + 3im + 2jm + 5km)
2
```
"""
jmag(q::Quaternion) = q.jm

"""
    kmag(q)
Return the k-imaginary part of the quaternion number `q`.
# Example
```jldoctest
julia> kmag(1 + 3im + 2jm + 5km)
5
```
"""
kmag(q::Quaternion) = q.km

real(::Type{Quaternion{T}}) where {T<:Real} = T

isreal(q::Quaternion) = iszero(q.im) && iszero(q.jm) && iszero(q.km)
isinteger(q::Quaternion) = isreal(q) && isinteger(q.re)
isfinite(q::Quaternion) = isfinite(q.re) && isfinite(q.im) && isfinite(q.jm) && isfinite(q.km)
isnan(q::Quaternion) = isnan(q.re) || isnan(q.im) || isnan(q.jm) || isnan(q.km)
isinf(q::Quaternion) = isinf(q.re) || isinf(q.im) || isinf(q.jm) || isinf(q.km)
iszero(q::Quaternion) = iszero(q.re) && iszero(q.im) && iszero(q.jm) && iszero(q.km)
iscomplex(q::Quaternion) = iszero(q.jm) && iszero(q.km)
isone(q::Quaternion) = isreal(q) && isone(q.re)

zero(::Type{Quaternion{T}}) where {T<:Real} = Quaternion{T}(zero(T), zero(T), zero(T), zero(T))

"""
    quat(r::Real [, i::Real, j::Real, k::Real])
    quat(c::Complex [, z::Complex])
Convert real numbers or arrays to quaternion. `i`, `j`, `k`, default to zero.
Convert a complex number to a quaternion, with `j` and `k` parts zero.
Convert a pair of complex numbers to a quaternion `c + z*j`.
Mostly just a more convienient way to call the constructor Quaternion(...).
"""
quat(a::Real, b::Real, c::Real, d::Real) = Quaternion(a, b, c, d)
quat(x::Real) = Quaternion(x)
quat(z::Complex) = Quaternion(z)
quat(q::Quaternion) = q
quat(z1::Complex, z2::Complex) = Quaternion(z1, z2)
quat(v::Vector{<:Real}) = Quaternion(v)

"""
    complex(q::Quaternion)
Returns the complex part of a quaternion, equivalent to `complex(q.re, q.im)`.
"""
complex(::Type{Quaternion{T}}) where {T<:Real} = Complex{T}
complex(q::Quaternion) = complex(q.re, q.im)

"""
    quat(T::Type)
Returns an appropriate type which can represent a value of type `T` as a quaternion.
Equivalent to `typeof(quat(zero(T)))`.
"""
quat(::Type{T}) where {T<:Real} = Quaternion{T}
quat(::Type{Complex{T}}) where {T<:Real} = Quaternion{T}
quat(::Type{Quaternion{T}}) where {T<:Real} = Quaternion{T}

vec(q::Quaternion) = vcat(q.re, q.im, q.jm, q.km)

function show(io::IO, q::Quaternion)
    r, i, j, k = vec(q)
    compact = get(io, :compact, false)
    show(io, r)
    for (x, xm) in [(i, "im"), (j, "jm"), (k, "km")]
        if signbit(x) && !isnan(x)
            x = -x
            print(io, compact ? "-" : " - ")
        else
            print(io, compact ? "+" : " + ")
        end
        show(io, x)
        if !(isa(x,Integer) && !isa(x,Bool) || isa(x,AbstractFloat) && isfinite(x))
            print(io, "*")
        end
        print(io, xm)
    end
end

function show(io::IO, q::Quaternion{Bool})
    if q == im
        print(io, "im")
    elseif q == jm
        print(io, "jm")
    elseif q == km
        print(io, "km")
    else
        print(io, "Quaternion($(q.re),$(q.im),$(q.jm),$(q.km))")
    end
end

function show(io::IO, ::MIME"text/html", q::Quaternion)
    r, i, j, k = vec(q)
    compact = get(io, :compact, false)
    show(io, r)
    for (x, xm) in [(i, "i"), (j, "j"), (k, "k")]
        if signbit(x) && !isnan(x)
            x = -x
            print(io, compact ? "-" : " - ")
        else
            print(io, compact ? "+" : " + ")
        end
        show(io, x)
        if !(isa(x,Integer) && !isa(x,Bool) || isa(x,AbstractFloat) && isfinite(x))
            print(io, "*")
        end
        print(io, "" * xm * "")
    end
end

function show(io::IO, ::MIME"text/html", q::Quaternion{Bool})
    if q == im
        print(io, "i")
    elseif q == jm
        print(io, "j")
    elseif q == km
        print(io, "k")
    else
        print(io, "Quaternion($(q.re),$(q.im),$(q.jm),$(q.km))")
    end
end

function read(s::IO, ::Type{Quaternion{T}}) where {T}
    r = read(s,T)
    i = read(s,T)
    j = read(s,T)
    k = read(s,T)
    Quaternion{T}(r,i,j,k)
end
function write(s::IO, q::Quaternion)
    write(s,q.re,q.im,q.jm,q.km)
end

==(q::Quaternion, w::Quaternion) = q.re == w.re && q.im == w.im && q.jm == w.jm && q.km == w.km
==(q::Quaternion, z::Complex) = iscomplex(q) && q.im == imag(z) && q.re == real(z)
==(q::Quaternion, x::Real) = isreal(q) && q.re == x

isequal(q::Quaternion, w::Quaternion) = isequal(q.re,w.re) & isequal(q.im,w.im) &
                                        isequal(q.jm,w.jm) & isequal(q.km,w.km)

#TODO: hash

"""
    conj(q::Quaternion)
Compute the complex conjugate of a quaternion `z`.
"""
conj(q::Quaternion) = Quaternion(q.re,-q.im,-q.jm,-q.km)

"""
    abs(q::Quaternion)
Compute the absolute value (euclidean norm) of a quaternion.
"""
abs(q::Quaternion)  = sqrt(abs2(q))

"""
    abs2(q::Quaternion)
Compute the squared euclidean norm of a quaternion.
Equivalent to `q * conj(q)` but computationally faster and returns a real.
"""
abs2(q::Quaternion) = q.re*q.re + q.im*q.im + q.jm*q.jm + q.km*q.km

inv(q::Quaternion)  = conj(q)/abs2(q)
inv(q::Quaternion{<:Integer}) = inv(float(q))

-(q::Quaternion) = Quaternion(-q.re, -q.im, -q.jm, -q.km)
+(q::Quaternion, w::Quaternion) = Quaternion(q.re + w.re, q.im + w.im,
                                             q.jm + w.jm, q.km + w.km)
-(q::Quaternion, w::Quaternion) = Quaternion(q.re - w.re, q.im - w.im,
                                             q.jm - w.jm, q.km - w.km)
*(q::Quaternion, w::Quaternion) = Quaternion(q.re*w.re - q.im*w.im - q.jm*w.jm - q.km*w.km,
                                             q.re*w.im + q.im*w.re + q.jm*w.km - q.km*w.jm,
                                             q.re*w.jm - q.im*w.km + q.jm*w.re + q.km*w.im,
                                             q.re*w.km + q.im*w.jm - q.jm*w.im + q.km*w.re)

+(x::Real, q::Quaternion) = Quaternion(x + q.re, q.im, q.jm, q.km)
+(q::Quaternion, x::Real) = Quaternion(q.re + x, q.im, q.jm, q.km)
+(z::Complex, q::Quaternion) = Quaternion(real(z) + q.re, imag(z) + q.im, q.jm, q.km)
+(q::Quaternion, z::Complex) = Quaternion(q.re + real(z), q.im + imag(z), q.jm, q.km)
function -(x::Real, q::Quaternion)
    # we don't want the default type for -(Bool)
    re = x - q.re
    Quaternion(re, - oftype(re, q.im), - oftype(re, q.jm), - oftype(re, q.km))
end
function -(z::Complex, q::Quaternion)
    # we don't want the default type for -(Bool)
    re = real(z) - q.re
    Quaternion(re, imag(z) - q.im, - oftype(re, q.jm), - oftype(re, q.km))
end
-(q::Quaternion, x::Real) = Quaternion(q.re - x, q.im, q.jm, q.km)
-(q::Quaternion, z::Complex) = Quaternion(q.re - real(z), q.im - imag(z), q.jm, q.km)
*(x::Real, q::Quaternion) = Quaternion(x * q.re, x * q.im, x * q.jm, x * q.km)
*(q::Quaternion, x::Real) = Quaternion(q.re * x, q.im * x, q.jm * x, q.km * x)

*(q::Quaternion, z::Complex) = Quaternion(q.re*z.re - q.im*z.im, q.re*z.im + q.im*z.re,
                                                      q.jm*z.re + q.km*z.im, q.km*z.re - q.jm*z.im)
*(z::Complex, q::Quaternion) = Quaternion(z.re*q.re - z.im*q.im, z.re*q.im + z.im*q.re,
                                                      z.re*q.jm - z.im*q.km, z.re*q.km + z.im*q.jm)

/(a::R, q::S) where {R<:Real,S<:Quaternion} = (T = promote_type(R,S); a*inv(T(q)))
/(q::Quaternion, x::Real) = Quaternion(q.re/x, q.im/x, q.jm/x, q.km/x)
/(q::Quaternion, w::Quaternion) = q * inv(w)

rand(r::AbstractRNG, ::SamplerType{Quaternion{T}}) where {T<:Real} =
    Quaternion(rand(r,T), rand(r,T), rand(r,T), rand(r,T))

"""
When the type argument is quaternion, the values are drawn
from the circularly symmetric quaternion normal distribution.
This is std normal in each component but with variance scaled by 1/4.
"""
randn(r::AbstractRNG, ::Type{Quaternion{T}}) where {T<:AbstractFloat} =
    Quaternion(T(0.5)*randn(r,T), T(0.5)*randn(r,T), T(0.5)*randn(r,T), T(0.5)*randn(r,T))

"""
    norm(q::Quaternion, p::Real=2)
Computes the `p` norm of the quaternion components of `q` as a vector.
For the default `p=2` euclidean norm, this is same as the quaternion `abs()`.
"""
norm(q::Quaternion{T}, p::Real=2) where {T<:Real} = norm(vec(q), p)

"""
    normalize(q::Quaternion, p::Real=2)
Computes the normalized unit quaternion. Optional choice of `p` norm.
"""
normalize(q::Quaternion{T}, p::Real=2) where {T<:Real} = Quaternion(normalize(vec(q), p))

exp(q::Quaternion{<:Integer}) = exp(float(q))
log(q::Quaternion{<:Integer}) = log(float(q))

### See Neil Dantam - Quaternion Computation - for exp and log
function exp(q::Quaternion{<:AbstractFloat})
    V = norm([q.im q.jm q.km])
    s = if V < 1e-8
        1 - V*V/6 + V*V*V*V/120
    else
        sin(V)/V
    end
    exp(q.re)*(cos(V) + quat(zero(q.re), q.im*s, q.jm*s, q.km*s))
end

function log(q::Quaternion{<:AbstractFloat})
    V = norm([q.im q.jm q.km])
    Q = abs(q)
    ϕ = atan(V, q.re)
    s = if V < 1e-8
        (1+ϕ*ϕ/6 +7*ϕ*ϕ*ϕ*ϕ/360)/Q
    else
        ϕ/V
    end
    quat(log(Q), q.im*s, q.jm*s, q.km*s)
end

### NOTE: sqrt(q^2) == q is NOT true for pure quaternion for this computation
sqrt(q::Quaternion) = exp(0.5*log(q))
^(q::Quaternion, w::Quaternion) = exp(w*log(q))

function round(q::Quaternion, rr::RoundingMode=RoundNearest, ri::RoundingMode=rr; kwargs...)
    Quaternion(round(real(q), rr; kwargs...),
               round(imag(q), ri; kwargs...),
               round(jmag(q), ri; kwargs...),
               round(kmag(q), ri; kwargs...))
end

float(q::Quaternion{<:AbstractFloat}) = q
float(q::Quaternion) = Quaternion(float(q.re), float(q.im), float(q.jm), float(q.km))

big(::Type{Quaternion{T}}) where {T<:Real} = Quaternion{big(T)}
big(q::Quaternion{T}) where {T<:Real} = Quaternion{big(T)}(q)

## Matrix representations of Quaternions

"""
    cmatrix(q::Quaternion)
    cmatrix(q::Matrix{Quaternion})
Returns the complex matrix representation of a quaternion.
If we write `q = x + y*j` then this is `[x y; -conj(y) conj(x)]`.
This works if `q` is a quaternion of matrix of quaternions.
For a matrix of size `(n,m)` this gives a complex matrix of size `(2n, 2m)`.
"""
function cmatrix(q::Quaternion)
    [complex( q.re, q.im) complex(q.jm,  q.km);
     complex(-q.jm, q.km) complex(q.re, -q.im)]
end

function cmatrix(Q::Matrix{Quaternion{T}}) where {T}
    [complex.( real.(Q),  imag.(Q)) complex.( jmag.(Q),  kmag.(Q));
     complex.(-jmag.(Q),  kmag.(Q)) complex.( real.(Q), -imag.(Q))]
end

"""
    qmatrix(c::Matrix{Complex})
Given a complex matrix representation of a quaternion of quaternion matrix,
returns the quaternion representation of the matrix of size `(n÷2, m÷2)`.
This gives an inverse of the `cmatrix` function.
This function does no error checking, make sure your matrix is valid beforehand.
"""
function qmatrix(C::Matrix{Complex{T}}) where {T}
    n, m = size(C)
    quat.(C[1:n÷2, 1:m÷2], C[1:n÷2, m÷2+1:m])
end

## Array operations on complex numbers ##

quat(A::AbstractArray{<:Quaternion}) = A


text:
TOL = 0.5e-8

function bisection(f::Function, a::Real, b::Real;
                   ε::Real = TOL, iₘₐₓ::Int = 100)

    fa = f(a)
    fb = f(b)
    if fa * fb >= 0
        throw(DomainError("The signs of the interval [f($(a)), f($(b)] must
                          differ."))
    end

    if a > b # enforces that the condition a ≦ b holds
        a, b = b, a
    end

    c = (a + b) / 2
    i = 0

    while (b - a) > 2ε && n < iₘₐₓ
        c = (a + b) / 2
        fc = f(c)
        if fc == 0
            break
        elseif fc * fa < 0
            b = c
            fb = fc
        else
            a = c
            fa = fc
        end
        i += 1
    end
    c, i
end


text:
#=
    Cheating model.

    The model is a more verbose implementation of the Cheating problem in
    https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter2_MorePyMC/Ch2_MorePyMC_PyMC3.ipynb
    """
    What we need is a better algorithm to ask students if they had cheated. Ideally the algorithm
    should encourage individuals to be honest while preserving privacy. The following proposed algorithm
    is a solution I greatly admire for its ingenuity and effectiveness:

        In the interview process for each student, the student flips a coin, hidden from the interviewer.
        The student agrees to answer honestly if the coin comes up heads. Otherwise, if the coin comes
        up tails, the student (secretly) flips the coin again, and answers "Yes, I did cheat" if the
        coin flip lands heads, and "No, I did not cheat", if the coin flip lands tails. This way, the
        interviewer does not know if a "Yes" was the result of a guilty plea, or a Heads on a second
        coin toss. Thus privacy is preserved and the researchers receive honest answers.
    """

    Compared to cheating.jl we add an array (cheated) which collected the students
    that we think cheated. Of course it's just a simulation since we cannot pinpoint
    the specific students that actually cheated.

    But the model works and give the same results as cheating.jl (and cheating2.jl), i.e. that
    it was about 22% of the student that cheated (the p variable).

        Summary Statistics
    parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

               p    0.2212    0.1192     0.0038    0.0063   307.7500    1.0006       28.8048
        obs_prop    0.3556    0.0435     0.0014    0.0028   377.7931    1.0154       35.3606
         p_prior    0.5006    0.2849     0.0090    0.0085   586.1921    1.0036       54.8664

   Cf cheating.jl and cheating2.jl for ports of the two PyMC3 models from the
      "Probabilistic Programming Hacker" book.

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function cheating_model0(N,X)
    p ~ Uniform(0,1) # freq cheating
    tail = 0
    head = 1

    cheated      ~ filldist(Bernoulli(p),N)
    truths       ~ filldist(Bernoulli(p),N)
    first_flips  ~ filldist(Bernoulli(0.5), N)
    second_flips ~ filldist(Bernoulli(0.5), N)

    # val = tzeros(N) # We must defined val to be used as LHS in @.

    # @. val = first_flips*truths + (1 - first_flips)*second_flips

    for i in 1:N
        # """
        # In the interview process for each student, the student flips a coin,
        # hidden from the interviewer.
        # """
        if first_flips[i] == head
            # The student agrees to answer honestly if the coin comes up heads.
            # """
            if cheated[i] == true
                truths[i] ~ Dirac(true)
            else
                truths[i] ~ Dirac(false)
            end
            second_flips[i] ~ Dirac(true)
        else
            # """
            # Otherwise, if the coin comes up tails, the student (secretly) flips the coin
            # again, and answers "Yes, I did cheat" if the coin flip lands heads,
            # and "No, I did not cheat", if the coin flip lands tails.
            # """
            if second_flips[i] == head
                truths[i] ~ Dirac(true)
            else
                truths[i] ~ Dirac(false)
            end
        end
    end

    obs_prop ~ Dirac(sum(truths) / N)
    X ~ Binomial(N, obs_prop )

    p_prior ~ Uniform(0,1)

end

N = 100
yes_response = 35
model = cheating_model0(N, yes_response)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5),  1_000)
# chns = sample(model, PG(5),  MCMCThreads(), 10_000, 4)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns[[:p,:obs_prop,:p_prior]])
display(plot(chns[[:p,:obs_prop,:p_prior]]))

# Another way:
println("mean(truths): $(mean(group(chns, :truths).value.data,dims=1) |> mean)")
println("mean(cheated): $(mean(group(chns, :cheated).value.data,dims=1) |> mean)")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#name: JuliaDup
#description: Duplicates a string in Julia
#language: julia
#tags: test, selenium
#input: string s
#output: string res
res = s^2

max_stars_repo_name: AndSviat/public

max_issues_repo_name: AndSviat/public

max_forks_repo_name: AndSviat/public


text:
## Farfield dynamics
#==========================================================================================#

"""
    farfield_velocity(r_i, r_j, Γ_j)

Compute the induced velocity at point ``r_i`` of the wake in the Trefftz plane due to the strength ``Γ_j`` at ``r_j``.
"""
farfield_velocity(r_i, r_j, Γ_j) = let r = r_i - r_j; Γ_j / 2π * SVector(1., 0., 0.) × r / norm(r)^2 end

"""
    farfield_influence_matrix(centers, normals, points)

Compute the aerodynamic influence coefficient matrix of the wake in the Trefftz plane given the center points, normal vectors, and the points of the wake panels.
"""
farfield_influence_matrix(centers, normals, points) = [ dot(farfield_velocity(r_i, r_j2, 1.) - farfield_velocity(r_i, r_j1, 1.), n_i) for (r_i, n_i) in zip(centers, normals), (r_j2, r_j1) in zip(points[2:end], points) ]

"""
    doublet_normal_derivatives(wake_points, Δφs, normals)

Compute the normal derivative strengths of the doublets given the wake points, net doublet strengths ``Δφ``s, and associated normal vectors.
"""
function doublet_normal_derivatives(wake_points, Δφs, normals)
    centers = @. (wake_points[1:end-1] + wake_points[2:end]) / 2
    AIC     = farfield_influence_matrix(centers, normals, wake_points)
    ∂φ_∂n   = AIC * Δφs
end

project_vector(vector, U) = vector - dot(U, vector) * U
normal(wake_proj_vector, U) = normalize(U × wake_proj_vector)
dihedral(wake_proj_vec)  = atan(wake_proj_vec[3], wake_proj_vec[2])

"""
    trefftz_plane_quantities(horseshoes :: AbstractArray{<: Horseshoe}, α, β)

Project `Horseshoe`s into the Trefftz plane aligned with the wind axes angles ``α,~β``, and compute normal vectors, projection angles and lengths.
"""
function trefftz_plane_quantities(horseshoes :: AbstractArray{<: Horseshoe}, α, β)
    # Reference velocity for broadcasting
    U_ref = (Ref ∘ SVector)(1, 0, 0)

    # Transform wake to wind axes
    wake_horsies = @views vec(horseshoes[end,:])                  # Trailing edge
    wake_points  = [ r1.(wake_horsies); [r2(wake_horsies[end])] ] # Trailing edge points
    wake_wind    = geometry_to_wind_axes.(wake_points, α, β)      # Wind axis transformation

    # Project trailing edge horseshoes' bound legs into Trefftz plane along wind axes
    wake_vectors   = @. wake_wind[2:end] - wake_wind[1:end-1] # Vectors
    wake_proj_vecs = @. project_vector(wake_vectors, U_ref)   # Project to wind axes

    # Normals, dihedral angles, and projection lengths
    wake_normals   = @. normal(wake_proj_vecs, U_ref)
    wake_dihedrals = @. dihedral(wake_proj_vecs)
    wake_lengths   = @. norm(wake_proj_vecs)

    wake_points, wake_normals, wake_dihedrals, wake_lengths
end

"""
    compute_farfield_forces(Δφs, Δs, ∂φ_∂n, θs, V, ρ)

Compute the aerodynamic forces in the Trefftz plane given cumulative doublet strengths ``Δφ``s, Trefftz panel lengths ``Δs``, doublet-normal directional derivatives ``∂φ/∂n``, Trefftz panel angles ``θ``s, the freestream speed ``V`` and density ``ρ``.
"""
function compute_farfield_forces(Δφs, Δs, ∂φ_∂n, θs, V, ρ)
    D_i = - 1/2 * ρ * sum(@. Δφs * Δs * ∂φ_∂n)
    Y   = - ρ * V * sum(@. Δφs * Δs * sin(θs))
    L   = ρ * V * sum(@. Δφs * Δs * cos(θs))

    SVector(D_i, Y, L)
end

"""
    farfield_forces(Γs, horseshoes, freestream, ρ)

Compute the aerodynamic forces in the Trefftz plane normal to the freestream given horseshoes, their associated strengths ``Γ``s, and a density ``ρ``.
"""
function farfield_forces(Γs, horseshoes, speed, α, β, ρ)
    # Get projections of horseshoes into Trefftz plane with the associated normals, dihedral angles and lengths
    wake_points, normals, dihedrals, Δs = trefftz_plane_quantities(horseshoes, α, β)

    # Compute directional derivatives of doublets in the normal direction
    Δφs   = vec(sum(Γs, dims = 1))
    ∂φ_∂n = doublet_normal_derivatives(wake_points, Δφs, normals)

    # Compute forces
    compute_farfield_forces(Δφs, Δs, ∂φ_∂n, dihedrals, speed, ρ)
end


text:
using Schistoxpkg
using Test
using Distributions
using Random
using JLD

egg_production_distribution = "NegBin"
N = 1000    #population size
time_step = 10.0
N_communities = 1
community_probs = [1.0]
community_contact_rate = [1.0]
density_dependent_fecundity = 0.0007
average_worm_lifespan = 5.7
max_age = 100
initial_worms = 0
initial_miracidia = 100000*N/1000
initial_miracidia_days = 3
init_env_cercariae = 100000*N/1000
worm_stages = 1
contact_rate = 0.1
max_fecundity = 0.34
max_fec_contact_rate_product =  0.034
contact_rate = max_fec_contact_rate_product / max_fecundity
female_factor = 1
male_factor = 1
age_contact_rates = [0.032,0.61, 1,0.06]
ages_for_contacts = [4,9,15,100]
contact_rate_by_age_array = fill(0.09, trunc(Int,(max_age+1)))
mda_adherence = 1
mda_access = 1
female_factor = 1
male_factor = 1
birth_rate = 28*time_step/(1000*365)
human_cercariae_prop = 1
predis_aggregation = 0.24
cercariae_survival = 1/2
miracidia_survival = 1/2
miracidia_maturity = 24
death_prob_by_age = [0.0656, 0.0093, 0.003, 0.0023, 0.0027, 0.0038, 0.0044, 0.0048, 0.0053,
                     0.0065, 0.0088, 0.0106, 0.0144, 0.021, 0.0333, 0.0529, 0.0851, 0.1366, 0.2183, 0.2998 , 0.3698, 1]

ages_for_death = [1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60,
                   65, 70, 75, 80, 85, 90, 95, 100, 110];
r = 0.03
vaccine_effectiveness  = 0.86
drug_effectiveness = 0.86
spec_ages = [7639, 7082, 6524, 5674, 4725, 4147, 3928, 3362,
              2636, 1970, 1468, 1166, 943, 718, 455, 244]
ages_per_index = 5
record_frequency = 1/24;
predis_aggregation = 0.8
scenario = "moderate adult";
use_kato_katz = 0
kato_katz_par = 0.87
heavy_burden_threshold = 16
rate_acquired_immunity = 0
M0 = 15
egg_sample_size  = 1
filename = "ee.jld"
human_larvae_maturity_time = 40

@testset "test_pars" begin
    pars = Parameters()
    @test pars.N == 1000
end

pars = Parameters(N, time_step, N_communities, community_probs, community_contact_rate,
        density_dependent_fecundity, average_worm_lifespan,
        max_age, initial_worms, initial_miracidia, initial_miracidia_days, init_env_cercariae,
        worm_stages, contact_rate, max_fec_contact_rate_product,  max_fecundity, age_contact_rates,
        ages_for_contacts, contact_rate_by_age_array, mda_adherence, mda_access,  female_factor, male_factor, miracidia_maturity,
        birth_rate, human_cercariae_prop, predis_aggregation, cercariae_survival, miracidia_survival,
        death_prob_by_age, ages_for_death, r, vaccine_effectiveness, drug_effectiveness,
        spec_ages, ages_per_index, record_frequency, use_kato_katz, kato_katz_par, heavy_burden_threshold, rate_acquired_immunity, M0,
        human_larvae_maturity_time, egg_sample_size, egg_production_distribution)
pars = make_age_contact_rate_array(pars, scenario, [],[]);

humans, miracidia, cercariae = create_population_specified_ages(pars)
@testset "test_human" begin
    h = Human()
    @test typeof(h) == Human
end

pars.cercariae_survival

@testset "create_contacts_settings" begin
    conts = create_contact_settings("low adult")
    @test conts[4] == 0.02
end

@testset "create_contacts_settings" begin
    conts = create_contact_settings("high adult")
    @test conts[4] == 0.12
end

@testset "create_contacts_settings" begin
    conts = create_contact_settings("moderate adult")
    @test conts[4] == 0.06
end

pars.max_age = 5

@test_throws ErrorException("max_age must be greater than 60") make_age_contact_rate_array(pars, "high adult", [], [])
pars.max_age = 100
pars = make_age_contact_rate_array(pars, "high adult", [], [])
@testset "contact_array" begin

    @test pars.contact_rate_by_age_array[1] == 0.01
end

humans, miracidia, cercariae = create_population_specified_ages(pars)

@testset "create_population" begin
    @test length(humans) == N
end

humans, miracidia, cercariae = create_population(pars)

@testset "create_population" begin
    @test length(humans) == N
end

humans = generate_ages_and_deaths(20000, humans, pars)
humans = update_contact_rate(humans,  pars)
@testset "update_contact_rate1" begin
    @test humans[1].age_contact_rate == pars.contact_rate_by_age_array[end]
end
humans[1].age = 0
humans = update_contact_rate(humans,  pars)

@testset "update_contact_rate2" begin
    @test humans[1].age_contact_rate == pars.contact_rate_by_age_array[1]
end

@testset "miracidia_death" begin
    @test isapprox(miracidia_death!(miracidia, pars), [pars.initial_miracidia, pars.initial_miracidia , round(pars.initial_miracidia  * pars.miracidia_survival)])
end

function calculate_time_step_specific_cerc_death(cercariae, pars, miracidia)
    updated_cercariae = cercariae
    time_step_specific_cerc_death = pars.cercariae_survival
    c = cercariae
    for i in 1:pars.time_step
        c = (c + miracidia[1]/pars.time_step) * pars.cercariae_survival
    end
    c1 = cercariae
    c1 = (c1 +  miracidia[1])
    time_step_specific_cerc_death = c/c1
    updated_cercariae = trunc(Int, round(cercariae * time_step_specific_cerc_death, digits= 0))
    return updated_cercariae
end
updated_cercariae = calculate_time_step_specific_cerc_death(cercariae, pars, miracidia)

@testset "cercariae_death" begin
    @test cercariae_death!(cercariae, miracidia, pars)== updated_cercariae

end

@testset "worm_pairs" begin
    @test isapprox(calculate_worm_pairs([[1,2],[2,2]],[[5,1],[2,0]]), [3,2])
end

humans[1].female_worms = [1,2]
humans[1].male_worms = [5,1]

female_worms = (p->p.female_worms).(humans)
male_worms = (p->p.male_worms).(humans)

@testset "worm_pairs" begin
    @test calculate_worm_pairs(female_worms, male_worms)[1] ==3
end

@testset "make_age_contact_rate_array(max_age,scenario)" begin
    @test pars.contact_rate_by_age_array[100] == 0.12
end

humans, miracidia, cercariae = create_population(pars)
old_cercariae = cercariae
old_length_miracidia = length(miracidia)

humans, cercariae, miracidia = cercariae_uptake!(humans, cercariae, miracidia, pars)
@testset "cerc_uptake" begin
    @test length(miracidia) <  old_length_miracidia
end

human_larvae = (p->p.larvae).(humans)

@testset "cerc_uptake" begin
    @test length(miracidia) <  old_length_miracidia
end

@testset "human_larvae" begin
    @test sum(sum.(human_larvae)) == 0
end

cercariae = 10000
humans, cercariae, miracidia  = cercariae_uptake_with_human_larvae!(humans, cercariae, miracidia, pars)

human_larvae = (p->p.larvae).(humans)
#
#
# @testset "human_larvae" begin
#     @test sum(sum.(human_larvae)) > 0
# end

human_larvae = (p->p.larvae).(humans)
x = sum(sum.(human_larvae))
humans = human_larvae_maturity(humans, pars)
human_larvae = (p->p.larvae).(humans)

@testset "human_larvae_maturity" begin
    @test sum(sum.(human_larvae)) == x
end

pars = make_age_contact_rate_array(pars, scenario, [4,9,15,pars.max_age], [10.1,2.1,4.2,6.2])
@testset "new_contact_rate" begin
    @test pars.contact_rate_by_age_array[1] == 10.1
end
start_worms = 10000
for h in humans
    h.female_worms = [start_worms]
    h.male_worms = [start_worms]
end
humans = worm_maturity!(humans, pars)
@testset "worm_maturity" begin
    @test humans[1].female_worms[1] < start_worms
end

#
# @testset "calc_worm_pairs" begin
#     @test isapprox(calculate_worm_pairs([[1,1],[2,3]], [[2,3],[5,3]]) , [2,5])
# end

pars.mda_adherence = 1
pars.mda_access = 1
pars.initial_worms = 50
Random.seed!(2525)
humans, miracidia, cercariae = create_population(pars)

humans = egg_production!(humans, pars)
# @testset "egg_prod" begin
#     @test mean((p->p.eggs).(humans)) == 5.402
# end
#
# @testset "mira_prod" begin
#     @test miracidia_production!(humans) == 81
# end
eggs = sum((p->p.eggs).(humans))
humans = death_of_human(humans)
@testset "death" begin
    length(humans) < pars.N
end

humans = administer_drug(humans, 1:length(humans), 1)
@testset "administer_drug" begin
    @test sum((p->p.eggs).(humans)) == eggs
end

@testset "administer_drug" begin
    @test sum(sum.((p->p.female_worms).(humans))) == 0
end

administer_vaccine(humans, 1:length(humans), 1, 3)

@testset "administer_vaccine" begin
    @test sum(sum.((p->p.vac_status).(humans))) == 3*length(humans)
end

Random.seed!(2525)
humans, miracidia, cercariae = create_population(pars)

humans = egg_production!(humans, pars)

mda(humans, 1, 0, 100, 1, [0,1])
@testset "mda" begin
    @test sum((p->p.eggs).(humans)) == 0
end

first_mda_time = 1
last_mda_time = 10
mda_info = create_mda(0, 0.75, 0, first_mda_time,
    last_mda_time, 1, [0,1], [0,1], [0,1], 1)
@testset "create_mda" begin
    @test length(mda_info) == (last_mda_time - first_mda_time)*3
end

mda_coverage, min_age_mda, max_age_mda, mda_effectiveness, next_mda_time, mda_gender = update_mda(mda_info, 0)
@testset "update_mda" begin
    @test [min_age_mda, max_age_mda] == [0,4]
end

@testset "update_mda" begin
 @test update_mda(mda_info, 1)[1] == 0.75
end

@testset "vac_decay" begin
    @test vac_decay!(humans,pars)[1].vac_status == -pars.time_step
end

Random.seed!(33)
@testset "Kato_katz" begin
    @test kato_katz(10, Gamma(0.87, 1/0.87)) == 1
end

Random.seed!(33)
humans, miracidia, cercariae = create_population(pars)
humans = egg_production!(humans, pars)

# @testset "count_eggs" begin
#     @test count_eggs(humans) == 5597
# end
#
#
# record = get_prevalences!(humans, time, pars)
# @testset "get_prevs" begin
#     @test isapprox(record.sac_burden, [93.94, 57.58, 1.01])
# end
#
save(filename, "humans", humans,  "miracidia", miracidia, "cercariae", cercariae, "pars", pars)
#
# humans_loaded,  miracidia_loaded, cercariae_loaded, pars_loaded = load_population_from_file(filename)
#
# @testset "load_pop" begin
#     @test humans_loaded[1].age == humans[1].age
# end

Random.seed!(33)

# @testset "update_env_to_equ" begin
#     @test update_env_to_equilibrium(1, humans, miracidia, cercariae, pars)[2][3] == 41
# end
#
#
# Random.seed!(33)
# @testset "update_env_to_equ_inc" begin
#     @test update_env_to_equilibrium_increasing(1, humans, miracidia, cercariae, pars)[2][3] == 14
# end
#
vaccine_info = []
mda_info = create_mda(0, .75, 0, 0, 2, 1, [0,1], [0,1], [0,1], pars.drug_effectiveness)
Random.seed!(33)

#
# @testset "update_env_const_pop" begin
#     @test isapprox(update_env_constant_population(10, humans,  miracidia, cercariae, pars, mda_info, vaccine_info)[2], [10160, 12146, 12620])
# end
#
#
#
#
# @testset "update_env_no_birth_death" begin
#     @test isapprox(update_env_no_births_deaths(10, humans,  miracidia, cercariae, pars, mda_info, vaccine_info)[2], [19110, 19696, 20028])
# end
#
#
#
#
#
# @testset "update_env_human_larvae" begin
# @test isapprox(update_env_to_equilibrium_human_larvae(5, humans, miracidia, cercariae, pars)[2],[20672, 20682, 20594])
# end
#
#
# @testset "update_env_human_larvae_constant_pop" begin
#     @test isapprox(update_env_constant_population_human_larvae(5, humans,  miracidia,
#                     cercariae, pars, mda_info, vaccine_info)[2],
#                         [23944, 24494, 25000])
# end
#
#
#
#
# @testset "update_env_no_births_deaths_human_larvae" begin
# @test isapprox(update_env_no_births_deaths_human_larvae(5, humans,  miracidia, cercariae, pars, mda_info, vaccine_info)[2],[27568, 28168, 28605])
# end
#
#
# @testset "update_env_no_births_deaths_increasing" begin
# @test isapprox(update_env_no_births_deaths_increasing(5, humans,  miracidia, cercariae, pars, mda_info, vaccine_info)[2],[870, 908, 888])
# end
#
# @testset "update_env_constant_population_increasing" begin
# @test isapprox(update_env_constant_population_increasing(5, humans,  miracidia, cercariae, pars, mda_info, vaccine_info)[2],[868, 888, 877])
# end

humans, miracidia, cercariae, record =
    update_env_constant_population(25, humans,  miracidia, cercariae, pars, mda_info, vaccine_info);

times = []
prev = []
sac_prev = []
high_burden = []
high_burden_sac =[]
adult_prev = []
high_adult_burden = []

times, prev, sac_prev, high_burden, high_burden_sac, adult_prev, high_adult_burden =
collect_prevs(times, prev, sac_prev, high_burden, high_burden_sac, adult_prev, high_adult_burden, record, 1)

times, prev, sac_prev, high_burden, high_burden_sac, adult_prev, high_adult_burden =
collect_prevs(times, prev, sac_prev, high_burden, high_burden_sac, adult_prev, high_adult_burden, record, 2)

#
# @testset "run_repeated_sims_no_population_change" begin
# @test run_repeated_sims_no_population_change(filename, 5, mda_info, vaccine_info, 1)[2][1][1] == 64.6
# end
#
# @testset "run_repeated_sims_no_population_change" begin
# @test run_repeated_sims_no_population_change_increasing(filename, 5, mda_info, vaccine_info, 1)[2][1][1] == 60.8
# end
#
# @testset "run_repeated_sims_no_population_change" begin
# @test run_repeated_sims_no_births_deaths(filename, 5, mda_info, vaccine_info, 1)[2][1][1] == 86.8
# end
#
# @testset "run_repeated_sims_no_births_deaths_human_larvae" begin
# @test run_repeated_sims_no_births_deaths_human_larvae(filename, 5, mda_info, vaccine_info, 1)[2][1][1] == 80.1
# end
#
#
# @testset "run_repeated_sims_no_births_deaths_increasing" begin
# @test run_repeated_sims_no_births_deaths_increasing(filename, 5, mda_info, vaccine_info, 1)[2][1][1] == 75.2
# end