text:
# SolarIrradiance.jl
# Copyright (c) 2021 Pierre Galloo-Beauvais, Pierre Haessig
# This code is subject to the terms of the MIT license (see LICENSE.txt)

module SolarIrradiance

export dayofyear, declination, EoT, hour_angle, day_bounds
export cosθavg, cosθZavg, Rb, global_radiation_tilt

include("GHIProcess.jl")
export GHIProcess

using Dates

"""
    declination(n)

Compute the solar declination angle δ of the day `n`, in radians.

δ is the angle between the Sun-Earth direction and the equatorial plane.
It varies between −23.45° on December 21-22 and +23.45° on June 22
and is zero at the equinoxes.
See https://www.pveducation.org/pvcdrom/properties-of-sunlight/declination-angle.

# Examples

Declination at the summer solstice in the northern hemisphere:
```jldoctest
julia> using Dates

julia> n = dayofyear(2021,06,22)
173

julia> δ = declination(n)
0.40924559967108437

julia> round(rad2deg(δ); digits=3)
23.448
```

Declination at March equinoxe:
```jldoctest
julia> n = dayofyear(2021,03,22)
81

julia> round(declination(n); digits=3)
-0.0
```

"""
declination(n) = deg2rad(23.45) * sin(2π*(284+n)/365)

"""
    EoT(n)

Equation of Time on day `n`, in hours.

The Equation of Time is the difference between the apparent (or true) solar time
and the mean solar time.
See https://www.pveducation.org/pvcdrom/properties-of-sunlight/solar-time

# Examples

Minimal value is about -14 minutes in mid February:
```jldoctest
julia> EoT(45)*60
-14.26116933144462
```

Maximal value is about +16 minutes in mid November:
```jldoctest
julia> EoT(306)*60
16.388668719846564
```
"""
function EoT(n)
    B= 2π*(n-1)/365
    return 3.82*(0.000075+0.001868*cos(B)-0.032077*sin(B)-0.014615*cos(2*B)-0.04089*sin(2*B))
end

"""
    hour_angle(n, tc, lon)

Compute the hour angle of the Sun ω (in rad) a civil hour `tc` on day `n`, at longitude `lon` (°).

See also the inverse function `time_hour_angle`.

L'angle horaire du Soleil donne un renseignement sur le temps solaire vrai, il vaut 0 quand le Soleil est à son point culminant de la journée, est négatif avant ce point et positif après ce point.

TODO: add missing time zone input
"""
hour_angle(n, tc, lon) = (tc + lon/15 + EoT(n) - 12)*deg2rad(15)

"""
    time_hour_angle(n, ω, lon)

Compute the civil hour `tc` of an hour angle of the Sun `ω` (in rad) on day `n`,
at longitude `lon` (°).

It is the inverse function of `tc -> hour_angle(n, tc, lon)`.

TODO: add missing time zone input
"""
time_hour_angle(n, ω, lon) = ω/deg2rad(15) - EoT(n) + 12 - lon/15

"""
    day_bounds(n, slope, azimuth, lat)

Compute the hour angle bounds when the sun is above a panel during day `n`.

Panel has a `slope` (in [0, 90°]) and `azimuth` (in [-90°, 90°])
and is located at latitude `lat` (°).

Return `ωu`, `ωd` (in rad) which are respectively the hour angles of
sunrise and sunset *over the panel*.
When `slope` is 0, this corresponds to regulat sunrise and sunset (over the ground).
and there is a dedicated method for this, see `day_bounds(n, lat)`.

Mathematically, with θ the angle of incidence of the Sun on the panel,
this means that cosθ=0 at hour angle `ωu` and `ωd` and cosθ ≥ 0 in between.
"""
function day_bounds(n, slope, azimuth, lat)
    beta=deg2rad(slope)
    gamma=deg2rad(azimuth)
    phi=deg2rad(lat)
    delta=declination(n)

    sinδ=sin(delta)
    cosδ=cos(delta)
    sinϕ=sin(phi)
    cosϕ=cos(phi)
    cosβ=cos(beta)
    sinβ=sin(beta)
    cosγ=cos(gamma)

    A = sinδ*sinϕ*cosβ - sinδ*cosϕ*sinβ*cosγ
    B = cosδ*cosϕ*cosβ + cosδ*sinϕ*sinβ*cosγ
    C = cosδ*sinβ*sin(gamma)

    D = C^2+B^2-A^2
    E = A-B
    if D >= 0
        ωu = 2*atan((-C+sqrt(D))/E)
        ωd = 2*atan((-C-sqrt(D))/E)
    else
        # quick fix to  avoid sqrt(negative number) error,
        # but there may be a more physical way to detect this situation
        ωu = 0.0
        ωd = 0.0
    end

    return ωu, ωd
end

"""
    day_bounds(n, lat)

Compute the hour angle bounds when the sun is above the horizon during day `n`
at latitude `lat` (°).

Return `ωu`, `ωd` (in rad) which are respectively the hour angles of sunrise and sunset.

Mathematically, with θz the solar zenith angle (angle of incidence of the Sun on the ground),
this means that cosθz=0 at hour angle `ωu` and `ωd` and cosθz ≥ 0 in between.
"""
function day_bounds(n, lat)
    phi=deg2rad(lat)

    # Sunset
    ωd = acos(-tan(phi)*tan(declination(n)))
    # Sunrise
    ωu = -ωd

    return ωu, ωd
end

"""
    cosθavg(n, tcini, dt, slope, azimuth, lat, lon)

Cosine of θ, the angle of incidence of the sun on a tilted panel.

Panel has orientation (`slope`, `azimuth` °) and is located at (`lon`, `lat` °).
Value is computed on day `n`, averaged on a time interval `dt`
starting at civil hour `tcini`.

See also: [`cosθZavg`](@ref) for the incidence on the _ground_.
"""
function cosθavg(n, tcini, dt, slope, azimuth, lat, lon)
    beta=deg2rad(slope)
    gamma=deg2rad(azimuth)
    phi=deg2rad(lat)

    ω1=hour_angle(n, tcini, lon)
    ω2=hour_angle(n, tcini+dt, lon)

    ωup, ωdown = day_bounds(n, slope, azimuth, lat)

    if ω1<ωup && ω2<ωup
        return 0.0
    elseif ω1>ωdown && ω2>ωdown
        return 0.0
    elseif ω1<ωup && ω2>ωup
        ω1 = ωup
        tdebut = time_hour_angle(n, ωup, lon)
        tfin = tcini+dt
        dt = tfin-tdebut # TODO: check that the change of `dt` is correct
    elseif ω1<ωdown && ω2>ωdown
        ω2 = ωdown
        tdebut = tcini
        tfin = time_hour_angle(n, ωdown, lon)
        dt = tfin-tdebut # TODO: check that the change of `dt` is correct
    end

    delta = declination(n)

    sinδ = sin(delta)
    cosδ = cos(delta)
    sinϕ = sin(phi)
    cosϕ = cos(phi)
    cosβ = cos(beta)
    sinβ = sin(beta)
    cosγ = cos(gamma)

    a = sinδ*sinϕ*cosβ - sinδ*cosϕ*sinβ*cosγ
    b = cosδ*cosϕ*cosβ + cosδ*sinϕ*sinβ*cosγ
    c = cosδ*sinβ*sin(gamma)
    f = ω2 - ω1
    g = sin(ω2) - sin(ω1)
    h = cos(ω2) - cos(ω1)
    x = (a*f + b*g - c*h)*12/(pi*dt)

    return x
end

"""
    cosθZavg(n, tcini, dt, lat, lon)

Cosine of ``θ_Z``, the solar zenith angle, i.e. the angle of incidence of the sun on the ground.

Location is (`lon`, `lat` °).
Value is computed on day `n`, averaged on a time interval `dt`
starting at civil hour `tcini`.

See also: [`cosθavg`](@ref) for the incidence on a _tilted panel_.
"""
function cosθZavg(n, tcini, dt, lat, lon)
    phi=deg2rad(lat)

    ω1=hour_angle(n, tcini, lon)
    ω2=hour_angle(n, tcini+dt, lon)

    ωup, ωdown = day_bounds(n, lat)

    if ω1<ωup && ω2<ωup
        return 0.0
    elseif ω1>ωdown && ω2>ωdown
        return 0.0
    elseif ω1<ωup && ω2>ωup
        ω1 = ωup
        tdebut = time_hour_angle(n, ωup, lon)
        tfin = tcini+dt
        dt = tfin-tdebut # TODO: check that the change of `dt` is correct
    elseif ω1<ωdown && ω2>ωdown
        ω2 = ωdown
        tdebut = tcini
        tfin = time_hour_angle(n, ωdown, lon)
        dt = tfin-tdebut # TODO: check that the change of `dt` is correct
    end

    delta = declination(n)

    d = cos(phi)*cos(delta)
    e = sin(phi)*sin(delta)
    f = ω2 - ω1
    g = sin(ω2) - sin(ω1)
    x = (d*g + e*f)*12/(pi*dt)

    return x
end

"""
    Rb(cosθ, cosθZ, Rbsat=5)

Rb = cosθ/cosθZ with special care of zeros and saturation at `Rbsat`
"""
function Rb(cosθ, cosθZ, Rbsat=5.)
    if cosθ==0. || cosθZ==0.
        return 0.0
    else
        x = cosθ/cosθZ
    end

    if x>Rbsat
        return Rbsat
    else
        return x
    end
end

"""
    Rb(n, tcini, dt, slope, azimuth, lat, lon, Rbsat=5)

Cette fonction calcule le terme de redressement de GHI direct pour un panneau PV incliné d'un angle `slope`, orienté d'un angle `azimuth`, situé à une latitude `lat` et une longitude `lon` (°). Elle le calcule pour un jour n donné, sur un intervalle donné de durée dt à partir d'un instant initial tcini.

⚠ Fonction avec encore quelques problèmes. Des pics apparaissent dans certaines conditions, une saturation a été mise en place afin que ces pics n'impactent pas les calculs utilisant cette fonction. Cette saturation, choisie empiriquement, empêche le terme de redressement de dépasser la valeur de 5.
"""
function Rb(n, tcini, dt, slope, azimuth, lat, lon, Rbsat=5.)
    cosθ = cosθavg(n, tcini, dt, slope, azimuth, lat, lon)
    cosθZ = cosθZavg(n, tcini, dt, lat, lon)

    return Rb(cosθ, cosθZ, Rbsat)
end

"""
    global_radiation_tilt(GHI::Real, n, tcini, dt, lat, lon, slope, azimuth, albedo;
                          diffuse_model=:EKD82, transpose_model=:HDKR)

Estimate the global irradiance on a tilted panel from horizontal data `GHI`.

Panel has a orientation (`slope`, `azimuth`) and is located at (`lon`, `lat`).
`GHI` data is for day `n`, averaged on a time interval `dt`,
starting at civil hour `tcini` (UTC).

Reflection from surrounding area is modeled by `albedo` in [0,1].

# Empirical models

Many steps of the irradiance transposition from a horizontal plane to a tilted plane
are deduced from geometry. However, two aspects of the estimation are _empirical_,
so that several models are available:
1. splitting the global horizontal irradiance between direct and diffuse components
2. transposing the diffuse irradiance on a tilted plane

## 1. Diffuse fraction estimation

`diffuse_model` sets the model used to estimate the _fraction of diffuse radiation_
from the clearness index:
- `:EKD82`: Erbs, Klein & Duffie (1982). Estimation of the diffuse radiation fraction for
  hourly, daily and monthly-average global radiation. Solar Energy 28 (4), 293–302.
- `:OH77`: Orgill & Hollands (1977). Correlation equation for hourly diffuse radiation
  on a horizontal surface. Sol. Energy 19, 357–359.
- `:CR79`: Collares-Pereira & Rabl (1979), The average distribution of solar radiation
  —correlations between diffuse and hemispherical and between daily and hourly insolation values.
  Solar Energy 22, 155.

## 2. Diffuse irradiance transposition

`transpose_model` sets the model to estimate the diffuse radiation on a tilted plane.
Models differ on their assumption on how the diffuse irradiance is spatially distributed
on the sky hemisphere.
See (Demain 2013) for an introduction to those models and Duffie & Beckman's book for HDKR.
Models are classified between isotropic and anisotropic.

Isotropic models:
- `:LJ62`: Liu & Jordan (1962). Daily insolation on surfaces tilted towards the equator. ASHRAE, 53:526-41.
- `:Ko86`: Korokanis (1986). On the choice of the angle of tilt for south facing solar collectors in the Athens basin area. Solar Energy, 36:217-25.
- `:Ba02`: Badescu (2002). 3D isotropic approximation for solar diffuse irradiance on tilted
surfaces. Renewable Energy, 26:221-3.

Anisotropic models:
- `:Wi82`: Willmot (1982). On the climatic optimization of the tilt and azimuth of
  flat-plate solar collectors. Solar Energy, 28:205-16.
- `:Bu77`: Bugler (1977). The determination of hourly insolation on an inclined plane
  using a diffuse irradiance model based on hourly measured global horizontal insolation.
  Solar Energy, 19:477-91.
- `:MI83`: Ma & Iqbal (1983) Statistical comparison of models for estimating solar radiation
  on inclined surfaces. Solar Energy 31 (3), 313-317.
- `:Ha79`:  Hay (1979). Study of shortwave radiation on non-horizontal surfaces.
  Canadian Climate Centre, Report No. 79-12, Downsview, Ontario.
- `:TC77`:  Temps & Coulson (1977). Solar radiation incident upon slopes of different orientation.
  Solar Energy, 19:179-84.
- `:Kl79`: Klucher (1979). Evaluation of models to predict insolation on tilted surfaces.
  Solar Energy, 23:111-114.
- `:HDKR`: HDKR model (as named in Duffie & Beckman's book) used by HOMER, proposed by
  Reindl, Beckman & Duffie (1990). Evaluation of Hourly Tilted Surface Radiation Models.
  Solar Energy, 45, 9-17.
"""
function global_radiation_tilt(GHI::Real, n, tcini, dt, lat, lon, slope, azimuth, albedo;
                               diffuse_model=:EKD82, transpose_model=:HDKR)
    delta = declination(n)
    phi = deg2rad(lat)
    beta = deg2rad(slope)
    gamma = deg2rad(azimuth)

    # Extraterrestrial normal irradiance
    Gon = 1367*(1+0.033*cos(deg2rad((360*n)/365)))

    # Averaged horizontal extraterrestrial irradiance
    ω1 = hour_angle(n, tcini, lon)
    ω2 = hour_angle(n, tcini+dt, lon)

    # TODO: use the cosθZavg function here instead:
    # Goavg = Gon*cosθZavg, I believe
    a = cos(phi)*cos(delta) * (sin(ω2) - sin(ω1))
    b = sin(phi)*sin(delta) * (ω2 - ω1)
    Goavg = Gon*(a+b)*12/(pi*dt)
    if Goavg < 0
        Goavg = 0.0
    end

    # Clearness index `kt`
    if Goavg > 0
        kt = GHI/Goavg
            if kt > 1
                kt = 1.0
            end
    else
        kt = 0.0
    end

    # Diffuse irradiance, estimated from `kt` and GHI using various models
    if diffuse_model == :EKD82 # Erbs et al. 1982
        if kt<=0.22
            GHIdiffus = (1-0.09*kt)*GHI
        elseif kt>0.8
            GHIdiffus = 0.165*GHI
        else
            GHIdiffus = (0.9511-0.1604*kt+4.388*kt^2-16.638*kt^3+12.336*kt^4)*GHI
        end

    elseif diffuse_model == :OH77 # Orgill & Hollands 1977
        if kt<=0.35
            GHIdiffus = (1-0.249*kt)*GHI
        elseif kt>0.75
            GHIdiffus = 0.177*GHI
        else
            GHIdiffus = (1.557-1.84*kt)*GHI
        end

    elseif diffuse_model == :CR79 # Colarres-Pereira & Rabl 1979
        if kt<=0.17
            GHIdiffus = 0.99*GHI
        elseif kt>=0.80
            GHIdiffus = 0.2*GHI
        elseif kt>0.17 && kt<0.75
            GHIdiffus = (1.188-2.272*kt+9.473*kt^2-21.865*kt^3+14.648*kt^4)*GHI
        else # [0.75, .80[ (not in original publication)
            GHIdiffus = (0.632-0.54*kt)*GHI
        end
    else
        error("unknown diffuse fraction estimation model $diffuse_model")
    end

    # Beam (direct) irradiance:
    GHIdirect=GHI-GHIdiffus

    # Direct (beam) radiation on tilted surface
    # averaged cos(θ) and cos(θz)
    cosθ = cosθavg(n, tcini, dt, slope, azimuth, lat, lon)
    cosθZ = cosθZavg(n, tcini, dt, lat, lon)
    RB=Rb(cosθ, cosθZ)
    Gdirect=RB*GHIdirect

    # Diffuse radiation on tilted surface
    if transpose_model == :LJ62 # Liu & Jordan (1962)
        Gdiffus=(1+cos(beta))*GHIdiffus/2

    elseif transpose_model == :Ko86 # Korokanis (1986)
        Gdiffus=(2+cos(beta))*GHIdiffus/3

    elseif transpose_model == :Ba02 # Badescu (2002)
        Gdiffus=(3+cos(2*beta))*GHIdiffus/4

    elseif transpose_model == :Wi82 # Willmot (1982)
        if cosθ <= 0
            AW = 0.0
            BW = 1.0115-0.20293*beta-0.080823*beta^2
        else
            AW = Gdirect*RB/(1367*cosθ)
            BW = (1.0115-0.20293*beta-0.080823*beta^2)*(1-(Gdirect/(1367*cosθ)))
        end
        Gdiffus = (AW+BW) * GHIdiffus

    elseif transpose_model == :Bu77 # Bugler (1977)
        AB = (1+cos(beta))/2
        if GHIdiffus <= 0
            BB = 0.0
        else
            BB = 0.05*(Gdirect/GHIdiffus)
        end
        if cosθZ <= 0
            CB = cosθ
        else
            CB = cosθ - (1+cos(beta))/(2*cosθZ)
        end
        Gdiffus = (AB + BB*CB) * GHIdiffus

    elseif transpose_model == :MI83 # Ma & Iqbal (1983)
        AM = (1-kt)*(1+cos(beta))/2
        BM = kt*RB
        Gdiffus = (AM+BM)*GHIdiffus

    elseif transpose_model == :Ha79 # Hay (1979)
        if Goavg <= 0
            AH = 1.0
            CH = 0.0
        else
            AH = 1 - (GHIdirect/Goavg)
            CH = RB * (GHIdirect/Goavg)
        end
        BH = (1+cos(beta))/2
        Gdiffus = (AH*BH+CH) * GHIdiffus

    elseif transpose_model == :TC77 # Temps & Coulson (1977)
        ATC = (1+cos(beta))/2
        BTC = 1+(cosθ^2)*(sqrt(1-cosθZ^2))^3
        CTC = 1+(sin(beta/2))^3
        Gdiffus = ATC*BTC*CTC*GHIdiffus

    elseif transpose_model == :Kl79 # Klucher (1979)
        AK = (1+cos(beta))/2
        if GHI <= 0
            FK = 1.0
        else
            FK = 1 - (GHIdiffus/GHI)^2
        end
        BK = 1 + FK*(cosθ^2)*(sqrt(1-cosθZ^2))^3
        CK = 1 + FK*(sin(beta/2))^3
        Gdiffus = AK*BK*CK * GHIdiffus

    elseif transpose_model == :HDKR # HDKR model from Reindl, Beckman & Duffie (1990)
        if Goavg > 0
            Ai = GHIdirect/Goavg
        else
            Ai = 0.0
        end
        if GHI > 0
            f = sqrt(GHIdirect/GHI)
        else
            f = 0.0
        end
        AHDKR = Ai*RB
        BHDKR = 1 - Ai
        CHDKR = (1+cos(beta))/2
        DHDKR = 1 + f*(sin(beta/2))^3
        Gdiffus = (AHDKR + BHDKR*CHDKR*DHDKR) * GHIdiffus

    else
        error("unknown diffuse irradiance transposition model $transpose_model")
    end

    # Reflected irradiance
    Greflechi = GHI * albedo * (1-cos(beta))/2

    # Total irradiance
    Gtot = Gdiffus + Gdirect + Greflechi

    # Saturate output
    if Gtot>2000.
        # TODO: transform fixed threshold into parameter
        Gtot = 2000.
    elseif Gtot<0
        Gtot = 0.
    end

    return Gtot
end

"""
    global_radiation_tilt(n, tcini, dt, lat, lon, fonctionGHI, diffuse_model, transpose_model, GHI_data, slope, azimuth, albedo)

Same as first `global_radiation_tilt(GHI::Real...)`, but

`GHI_data` is the Vector of hourly GHI of a given year, obtained from function `GHIannee`.

Elle prend également en argument une chaîne de caractères `fonctionGHI`, qui permet de choisir l'interpolation des données du GHI. Cette chaîne de caractère peut donc être : "forwardfill", "backwardfill", "centeredfill" et "interpolation".
"""
function global_radiation_tilt(GHI_data::Vector, fonctionGHI, n, tcini, dt, lat, lon, slope, azimuth, albedo;
                               diffuse_model=:EKD82, transpose_model=:HDKR)
    GHI = GHIProcess.GHIavg(GHI_data, n, tcini, dt, fonctionGHI)
    Gtot = global_radiation_tilt(GHI, n, tcini, dt, lat, lon, slope, azimuth, albedo;
                                 diffuse_model=:EKD82, transpose_model=:HDKR)
    return Gtot
end

end # module


text:
using HCIToolbox
using Documenter

setup = quote
    using HCIToolbox
    using Random
    Random.seed!(8799)
end

DocMeta.setdocmeta!(HCIToolbox, :DocTestSetup, setup; recursive = true)

makedocs(;
    modules = [HCIToolbox],
    authors = "Miles Lucas <mdlucas@hawaii.edu>",
    repo = "https://github.com/juliahci/HCIToolbox.jl/blob/{commit}{path}#L{line}",
    sitename = "HCIToolbox.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://juliahci.github.io/HCIToolbox.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
        "Processing" => "processing.md",
        "Spectral Processing" => "sdi.md",
        "Geometries" => "geometry.md",
        "Signal Injection" => "inject.md",
        "Utilites" => "utils.md",
        "Index" => "api.md"
    ],
)

deploydocs(;
    repo = "github.com/JuliaHCI/HCIToolbox.jl",
    push_preview = true
)

ext: jl


text:
@testset "Homotopies" begin
    @testset "StraightLineHomotopy" begin
        F = SPSystem(equations(katsura(5)))
        G = SPSystem(equations(cyclic(6)))
        H = StraightLineHomotopy(F, G)
        @test H isa AbstractHomotopy
        @test size(H) == (6, 6)
        @test HC.gamma(H) isa ComplexF64
        HomotopyContinuation.homotopy_interface_test(H)
    end

    @testset "HomotopyWithCache" begin
        x = rand(Complex{Float64}, 6)
        t = rand()
        F = SPSystem(equations(katsura(5)))
        G = SPSystem(equations(cyclic(6)))
        H = HomotopyWithCache(StraightLineHomotopy(F, G), x, t)
        @test H isa AbstractHomotopy
        @test size(H) == (6, 6)

        m, n = size(H)
        u = zeros(Complex{Float64}, m)
        U = zeros(Complex{Float64}, m, n)

        evaluate!(u, H, x, t)
        @test evaluate(H, x, t) ≈ u

        dt!(u, H, x, t)
        @test dt(H, x, t) ≈ u

        jacobian!(U, H, x, t)
        @test jacobian(H, x, t) ≈ U

        evaluate_and_jacobian!(u, U, H, x, t)
        @test all(evaluate_and_jacobian(H, x, t) .≈ (u, U))
        @test all((u, U) .≈ (evaluate(H, x, t), jacobian(H, x, t)))

        jacobian_and_dt!(U, u, H, x, t)
        @test all(jacobian_and_dt(H, x, t) .≈ (U, u))
        @test all((U, u) .≈ (jacobian(H, x, t), dt(H, x, t)))
    end

    @testset "ParameterHomotopy" begin
        F = equations(katsura(5))
        vars = variables(F)
        H = ParameterHomotopy(F, vars[5:6], p₁ = rand(2), p₀ = rand(2))
        @test H isa AbstractHomotopy
        @test size(H) == (6, 4)
        @test size(H, 1) == 6
        @test size(H, 2) == 4
        @test length(H) == 6
        @test HC.nvariables(H) == 4

        HomotopyContinuation.homotopy_interface_test(H)

        H = ParameterHomotopy(
            F,
            vars[5:6],
            p₁ = rand(2),
            p₀ = rand(2),
            γ₁ = randn(ComplexF64),
            γ₀ = randn(ComplexF64),
        )
        @test H isa ParameterHomotopy
        HomotopyContinuation.homotopy_interface_test(H)
    end

    @testset "FixedPointHomotopy" begin
        F = SPSystem(equations(katsura(5)))
        H = FixedPointHomotopy(F, rand(ComplexF64, 6))
        @test H isa AbstractHomotopy
        @test size(H) == (6, 6)

        HomotopyContinuation.homotopy_interface_test(H)
    end

    @testset "PatchedHomotopy" begin
        F = SPSystem(equations(katsura(5)))
        G = SPSystem(equations(cyclic(6)))
        x = ProjectiveVectors.embed(rand(Complex{Float64}, 5))
        H = PatchedHomotopy(StraightLineHomotopy(F, G), OrthogonalPatch(), x)
        @test H isa AbstractHomotopy
        @test size(H) == (7, 6)

        HomotopyContinuation.homotopy_interface_test(H, x)
    end

    @testset "Coefficient Homotopy" begin
        E = [[2 1 0; 0 0 0], [1 0; 1 0]]
        start = [[1.0 + 0im, -3.0, 2.0], [2.0 + 0im, -2.0]]
        target = [[2.0 + 0im, -2.0, 5.0], [3.0 + 0im, -1.0]]
        H = CoefficientHomotopy(E, start, target)
        x = randn(ComplexF64, 2)
        HomotopyContinuation.homotopy_interface_test(H, x)
    end

    @testset "ConstantHomotopy" begin
        F = SPSystem(equations(katsura(5)))
        H = ConstantHomotopy(F)
        @test H isa AbstractHomotopy
        @test size(H) == (6, 6)

        HomotopyContinuation.homotopy_interface_test(H)
    end
end


text:
struct PriorAM
  K::Int
  beta::Bool
  gamma::Beta
  eta::Dirichlet
  mu::Normal
  sigma::LogNormal
  nu::LogNormal
  phi::Normal
end

function PriorAM(K, beta; gamma=Beta(1,1), eta=nothing,
                 mu=Normal(0, 3), sigma=LogNormal(0, .5),
                 nu=LogNormal(1.6, 0.4), phi=Normal(-1, 10), data=nothing)
  eta == nothing && (eta = Dirichlet(K, 1/K))
  data == nothing || (mu = compute_prior_mu(data))
  return PriorAM(K, beta, gamma, eta, mu, sigma, nu, phi)
end

function logabsjacobian(v::AbstractVector{<:Real}, prior::PriorAM)
  K = prior.K

  out = MCMC.logabsjacobian_logitx(gammaCfromvec(K, v))
  out += MCMC.logabsjacobian_logitx(gammaTfromvec(K, v))

  out += MCMC.simplex_logabsdet(etaCfromvec(K, v))
  out += MCMC.simplex_logabsdet(etaTfromvec(K, v))
  out += sum(MCMC.logabsjacobian_logx.(sigmafromvec(K, v)))
  out += sum(MCMC.logabsjacobian_logx.(nufromvec(K, v)))

  return out
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
using DataStructures: DisjointSets, union!, in_same_set, num_groups

"Compute array of edges (Nodes in columns), sorted by length"
function sorted_edges(nodes::Vector{Node})
    nnodes = length(nodes)
    edges = []
    for i in 1:nnodes
        for j in i+1:nnodes
            push!(edges, nodes[[i,j]])
        end
    end
    sortslices(hcat(edges...), dims=2, by=edge -> begin
             u, v = edge
             (u.x - v.x)^2 + (u.y - v.y)^2
             end)
end

"""
Compute a minimum spanning tree using Kruskal's algorithm.

Returns a vector of arcs with integer indices for nodes.
"""
function min_spanning_tree(nodes::Vector{Node})
    # MAYBE, use a Delaunay triangulation instead?
    n = length(nodes)

    edges = sorted_edges(nodes)
    @assert size(edges, 1) == 2 # column-wise

    nodeidx = Dict(nodes[i] => i for i in 1:n)

    tree_edges = Arc[]
    sizehint!(tree_edges, n-1)

    sets = DisjointSets{Node}(nodes)
    for e in 1:size(edges, 2)
        u, v = edges[:, e]
        if !in_same_set(sets, u, v)
            union!(sets, u, v)
            push!(tree_edges, Arc(nodeidx[u], nodeidx[v]))
        end
        num_groups(sets) == 1 && break # spanning tree
    end

    tree_edges
end

"Build a Topology from an Instance's nodes, using MST"
function topology_from_mst(nodes::Vector{Node})
    Topology(nodes, min_spanning_tree(nodes))
end
topology_from_mst(instance::Instance) = topology_from_mst(instance.nodes)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualNetworkGatewayPropertiesFormat <: SwaggerModel
    ipConfigurations::Any # spec type: Union{ Nothing, Vector{VirtualNetworkGatewayIPConfiguration} } # spec name: ipConfigurations
    gatewayType::Any # spec type: Union{ Nothing, String } # spec name: gatewayType
    vpnType::Any # spec type: Union{ Nothing, String } # spec name: vpnType
    vpnGatewayGeneration::Any # spec type: Union{ Nothing, String } # spec name: vpnGatewayGeneration
    enableBgp::Any # spec type: Union{ Nothing, Bool } # spec name: enableBgp
    enablePrivateIpAddress::Any # spec type: Union{ Nothing, Bool } # spec name: enablePrivateIpAddress
    activeActive::Any # spec type: Union{ Nothing, Bool } # spec name: activeActive
    gatewayDefaultSite::Any # spec type: Union{ Nothing, SubResource } # spec name: gatewayDefaultSite
    sku::Any # spec type: Union{ Nothing, VirtualNetworkGatewaySku } # spec name: sku
    vpnClientConfiguration::Any # spec type: Union{ Nothing, VpnClientConfiguration } # spec name: vpnClientConfiguration
    bgpSettings::Any # spec type: Union{ Nothing, BgpSettings } # spec name: bgpSettings
    customRoutes::Any # spec type: Union{ Nothing, AddressSpace } # spec name: customRoutes
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState
    enableDnsForwarding::Any # spec type: Union{ Nothing, Bool } # spec name: enableDnsForwarding
    inboundDnsForwardingEndpoint::Any # spec type: Union{ Nothing, String } # spec name: inboundDnsForwardingEndpoint

    function VirtualNetworkGatewayPropertiesFormat(;ipConfigurations=nothing, gatewayType=nothing, vpnType=nothing, vpnGatewayGeneration=nothing, enableBgp=nothing, enablePrivateIpAddress=nothing, activeActive=nothing, gatewayDefaultSite=nothing, sku=nothing, vpnClientConfiguration=nothing, bgpSettings=nothing, customRoutes=nothing, resourceGuid=nothing, provisioningState=nothing, enableDnsForwarding=nothing, inboundDnsForwardingEndpoint=nothing)
        o = new()
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("ipConfigurations"), ipConfigurations)
        setfield!(o, Symbol("ipConfigurations"), ipConfigurations)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("gatewayType"), gatewayType)
        setfield!(o, Symbol("gatewayType"), gatewayType)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("vpnType"), vpnType)
        setfield!(o, Symbol("vpnType"), vpnType)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("vpnGatewayGeneration"), vpnGatewayGeneration)
        setfield!(o, Symbol("vpnGatewayGeneration"), vpnGatewayGeneration)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("enableBgp"), enableBgp)
        setfield!(o, Symbol("enableBgp"), enableBgp)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("enablePrivateIpAddress"), enablePrivateIpAddress)
        setfield!(o, Symbol("enablePrivateIpAddress"), enablePrivateIpAddress)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("activeActive"), activeActive)
        setfield!(o, Symbol("activeActive"), activeActive)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("gatewayDefaultSite"), gatewayDefaultSite)
        setfield!(o, Symbol("gatewayDefaultSite"), gatewayDefaultSite)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("sku"), sku)
        setfield!(o, Symbol("sku"), sku)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("vpnClientConfiguration"), vpnClientConfiguration)
        setfield!(o, Symbol("vpnClientConfiguration"), vpnClientConfiguration)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("bgpSettings"), bgpSettings)
        setfield!(o, Symbol("bgpSettings"), bgpSettings)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("customRoutes"), customRoutes)
        setfield!(o, Symbol("customRoutes"), customRoutes)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("enableDnsForwarding"), enableDnsForwarding)
        setfield!(o, Symbol("enableDnsForwarding"), enableDnsForwarding)
        validate_property(VirtualNetworkGatewayPropertiesFormat, Symbol("inboundDnsForwardingEndpoint"), inboundDnsForwardingEndpoint)
        setfield!(o, Symbol("inboundDnsForwardingEndpoint"), inboundDnsForwardingEndpoint)
        o
    end
end # type VirtualNetworkGatewayPropertiesFormat

const _property_map_VirtualNetworkGatewayPropertiesFormat = Dict{Symbol,Symbol}(Symbol("ipConfigurations")=>Symbol("ipConfigurations"), Symbol("gatewayType")=>Symbol("gatewayType"), Symbol("vpnType")=>Symbol("vpnType"), Symbol("vpnGatewayGeneration")=>Symbol("vpnGatewayGeneration"), Symbol("enableBgp")=>Symbol("enableBgp"), Symbol("enablePrivateIpAddress")=>Symbol("enablePrivateIpAddress"), Symbol("activeActive")=>Symbol("activeActive"), Symbol("gatewayDefaultSite")=>Symbol("gatewayDefaultSite"), Symbol("sku")=>Symbol("sku"), Symbol("vpnClientConfiguration")=>Symbol("vpnClientConfiguration"), Symbol("bgpSettings")=>Symbol("bgpSettings"), Symbol("customRoutes")=>Symbol("customRoutes"), Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("enableDnsForwarding")=>Symbol("enableDnsForwarding"), Symbol("inboundDnsForwardingEndpoint")=>Symbol("inboundDnsForwardingEndpoint"))
const _property_types_VirtualNetworkGatewayPropertiesFormat = Dict{Symbol,String}(Symbol("ipConfigurations")=>"Vector{VirtualNetworkGatewayIPConfiguration}", Symbol("gatewayType")=>"String", Symbol("vpnType")=>"String", Symbol("vpnGatewayGeneration")=>"String", Symbol("enableBgp")=>"Bool", Symbol("enablePrivateIpAddress")=>"Bool", Symbol("activeActive")=>"Bool", Symbol("gatewayDefaultSite")=>"SubResource", Symbol("sku")=>"VirtualNetworkGatewaySku", Symbol("vpnClientConfiguration")=>"VpnClientConfiguration", Symbol("bgpSettings")=>"BgpSettings", Symbol("customRoutes")=>"AddressSpace", Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"ProvisioningState", Symbol("enableDnsForwarding")=>"Bool", Symbol("inboundDnsForwardingEndpoint")=>"String")
Base.propertynames(::Type{ VirtualNetworkGatewayPropertiesFormat }) = collect(keys(_property_map_VirtualNetworkGatewayPropertiesFormat))
Swagger.property_type(::Type{ VirtualNetworkGatewayPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkGatewayPropertiesFormat[name]))}
Swagger.field_name(::Type{ VirtualNetworkGatewayPropertiesFormat }, property_name::Symbol) =  _property_map_VirtualNetworkGatewayPropertiesFormat[property_name]

const _allowed_VirtualNetworkGatewayPropertiesFormat_gatewayType = ["Vpn", "ExpressRoute"]

const _allowed_VirtualNetworkGatewayPropertiesFormat_vpnType = ["PolicyBased", "RouteBased"]

const _allowed_VirtualNetworkGatewayPropertiesFormat_vpnGatewayGeneration = ["None", "Generation1", "Generation2"]

function check_required(o::VirtualNetworkGatewayPropertiesFormat)
    true
end

function validate_property(::Type{ VirtualNetworkGatewayPropertiesFormat }, name::Symbol, val)
    if name === Symbol("gatewayType")
        Swagger.validate_param(name, "VirtualNetworkGatewayPropertiesFormat", :enum, val, _allowed_VirtualNetworkGatewayPropertiesFormat_gatewayType)
    end
    if name === Symbol("vpnType")
        Swagger.validate_param(name, "VirtualNetworkGatewayPropertiesFormat", :enum, val, _allowed_VirtualNetworkGatewayPropertiesFormat_vpnType)
    end
    if name === Symbol("vpnGatewayGeneration")
        Swagger.validate_param(name, "VirtualNetworkGatewayPropertiesFormat", :enum, val, _allowed_VirtualNetworkGatewayPropertiesFormat_vpnGatewayGeneration)
    end
end


text: include("test_SolveOneNonlinearEquation.jl")


text:
#=
Julia program to simulate spreading of COVID-19
Julia translation of https://github.com/MuAuan/collective_particles
=#

"""
Module for global parameters and variables
"""
module ParamVar
    struct Parameters
        num_particles::Int64  # Number of particles
        max_iteration::Int64  # Total number of time iteration
        #
        x_range::Float64  # Size of computational domain [0, x_range]
        y_range::Float64
        #
        vel_mean::Float64  # Average of velocity component
        vel_σ::Float64  # Standard deviation
        vel_flc::Float64  # Range of fluctuation
        #
        ratio_infection_init::Float64  # Initial percentage of infected particles
        recovery_time::Int64  # Recovery time
        infection_chance::Float64  # Chance of infection per unit time
        #
        radius_infection::Float64  # Radius of infection-zone
    end

    struct Flags
        flag_multiple_infection::Bool  # Multiple infection
        flag_infected_isolation::Bool  # Isolation (stop) of infected particle
    end

    mutable struct Particle
        pos_x::Float64  # Position x-component
        pos_y::Float64
        vel_x::Float64  # Velocity x-component
        vel_y::Float64
        status::AbstractString  # "not_infected" "infected" or "recovered"
        t_ifcn::Int64  # Infected time
        past_ifcn::Int64  # Past infection history
        # Constructor
        Particle() = new()
    end

    mutable struct NumSnapshot
        not_infected::Int64  # Number of never-infected particles
        infected::Int64  # Number of infected particles
        recovered::Int64  # Number of recovered particles
    end
end

"""
Module for simulation
"""
module TimeMarch
using Distributions
    """
    Compute initial condition
    """
    function set_initial_condition(param, ptcl)
        num_infected_init = Int64(param.ratio_infection_init * param.num_particles)
        for itr_ptcl = 1:param.num_particles
            ptcl[itr_ptcl].pos_x = rand(Uniform(0.0, param.x_range))  # Uniform distribution
            ptcl[itr_ptcl].pos_y = rand(Uniform(0.0, param.y_range))
            ptcl[itr_ptcl].vel_x = rand(Normal(param.vel_mean, param.vel_σ))  # Gaussian distribution
            ptcl[itr_ptcl].vel_y = rand(Normal(param.vel_mean, param.vel_σ))
            ptcl[itr_ptcl].t_ifcn = 0
            ptcl[itr_ptcl].past_ifcn = 0
            ptcl[itr_ptcl].status = "not_infected"  # Initially not infected
            # Some particles are initially infected
            if itr_ptcl <= num_infected_init
                ptcl[itr_ptcl].status = "infected"
                ptcl[itr_ptcl].past_ifcn = 1
            end
        end
        # ptcl.pos_x .= rand(Uniform(0.0, param.x_range))  # ERROR: LoadError: type Array has no field pos_x
    end

    """
    Compute (squared) relative distance of two particles
    """
    function compute_relative_distance(x1, y1, x2, y2)
        return (x1-x2)^2 + (y1-y2)^2
    end

    """
    Check whether the particle will be infected or not
    """
    function compute_is_infected(
        infection_chance, flag_multiple_infection,
        status, past_ifcn
    )
        isinfected = false

        # Infection chance & past infection criteria
        if rand(Uniform(0.0, 1.0)) <= infection_chance && status == "not_infected"
            isinfected = true
        end

        if flag_multiple_infection  # Multiple infection is allowed
            if status == "recovered"  # Infected in the past
                # Infection chance criteria is modified by past number of infection
                if rand(Uniform(0.0, 1.0)) <= infection_chance * (0.5^past_ifcn)
                    isinfected = true
                end
            end
        end

        return isinfected
    end

    """
    Compute updated velocity components
    """
    function compute_new_velocity(
        flag_infected_isolation, vel_flc, status,
        vx, vy
    )

        vx_flc = vel_flc * rand(Uniform(-1.0, 1.0))
        vy_flc = vel_flc * rand(Uniform(-1.0, 1.0))

        vx_new = vx + vx_flc
        vy_new = vy + vy_flc

        if flag_infected_isolation && status == "infected"
            vx_new, vy_new = 0.0, 0.0
        end

        return vx_new, vy_new
    end

    """
    Compute new coordinates of particles ensuring periodic boundary condition
    point ∈ [0, range]
    """
    function ensure_periodic_bc(point, range)
        if point > range
            new_point = point - range
        elseif point < 0.0
            new_point = point + range
        else
            new_point = point
        end

        return new_point
    end

    """
    Update particle properties
    """
    function update_particles(param, flag, ptcl)
        x_ifcn = Array{Float64}(undef, 0)
        y_ifcn = Array{Float64}(undef, 0)

        # Update properties of infected particles
        for itr_ptcl = 1:param.num_particles
            # Recovery from infection
            s = ptcl[itr_ptcl].status
            t = ptcl[itr_ptcl].t_ifcn
            if s == "infected" && t >= param.recovery_time
                ptcl[itr_ptcl].status = "recovered"  # Hold infection history
                ptcl[itr_ptcl].past_ifcn += 1  # Hold infection history
            end

            # Store position of infected particles
            if s == "infected"
                append!(x_ifcn, ptcl[itr_ptcl].pos_x)
                append!(y_ifcn, ptcl[itr_ptcl].pos_y)
                ptcl[itr_ptcl].t_ifcn += 1
            end
        end

        # Update properties of not_infected and recovered particles
        for itr_ptcl = 1:param.num_particles
            x = ptcl[itr_ptcl].pos_x
            y = ptcl[itr_ptcl].pos_y
            s = ptcl[itr_ptcl].status
            t = ptcl[itr_ptcl].t_ifcn
            p = ptcl[itr_ptcl].past_ifcn

            # Loop of infected particles
            for itr_ifcn = 1:length(x_ifcn)
                r2 = compute_relative_distance(x, y, x_ifcn[itr_ifcn], y_ifcn[itr_ifcn])
                if r2 < param.radius_infection^2  # Relative distance criteria
                    isinfected = compute_is_infected(
                        param.infection_chance, flag.flag_multiple_infection,
                        s, p)
                    if isinfected
                        s = "infected"
                        t = 0  # Time since infection
                        p += 1  # Past number of infection
                    end
                end
            end

            ptcl[itr_ptcl].status = s
            ptcl[itr_ptcl].t_ifcn = t
            ptcl[itr_ptcl].past_ifcn = p
        end

        # Update position and velocity of all particles
        for itr_ptcl = 1:param.num_particles
            x = ptcl[itr_ptcl].pos_x
            y = ptcl[itr_ptcl].pos_y
            vx = ptcl[itr_ptcl].vel_x
            vy = ptcl[itr_ptcl].vel_y
            s = ptcl[itr_ptcl].status

            # Update velocity
            vx_new, vy_new = compute_new_velocity(
                flag.flag_infected_isolation, param.vel_flc, s,
                vx, vy)

            # Update position
            x_new = x + vx_new
            y_new = y + vy_new

            # Ensure periodic boundary condition
            x_new = ensure_periodic_bc(x_new, param.x_range)
            y_new = ensure_periodic_bc(y_new, param.y_range)

            ptcl[itr_ptcl].pos_x = x_new
            ptcl[itr_ptcl].pos_y = y_new
            ptcl[itr_ptcl].vel_x = vx_new
            ptcl[itr_ptcl].vel_y = vy_new
        end
    end

    """
    Count number of
    - not infected particles
    - infected particles
    - recovered particles
    """
    function count_status(param, ptcl)
        num_not_infected, num_infected, num_recovered = 0, 0, 0

        for itr_ptcl = 1:param.num_particles
            # Count number of particles according to their status
            if ptcl[itr_ptcl].status == "not_infected"  # not infected
                num_not_infected += 1
            elseif ptcl[itr_ptcl].status == "infected"  # infected
                num_infected += 1
            elseif ptcl[itr_ptcl].status == "recovered"  # recovered
                num_recovered += 1
            else
                throw(DomainError(ptcl[itr_ptcl].status, "status must be 'not_infected', 'infected' or 'recovered'"))
            end
        end

        return num_not_infected, num_infected, num_recovered
    end
end

"""
Module for plot
"""
module Output
    using Printf
    using Base.Filesystem
    using Plots
    using StatsPlots
    font = Plots.font("Times New Roman", 20)

    """
    Stdout simulation condition
    """
    function stdout_condition(param, flag)
        println("#========== Coronavirus spreading simulation ====================")
        println(@sprintf "num_particles = %i,  max_iteration = %i" param.num_particles param.max_iteration)
        println(@sprintf "(x, y) = [0.00, %.2f] × [0.00, %.2f]" param.x_range param.y_range)
        println(@sprintf "vel_mean = %.3f,  vel_σ = %.3f,  vel_flc = %.3f" param.vel_mean param.vel_σ param.vel_flc)
        println(@sprintf "ratio_infection_init = %.3f" param.ratio_infection_init)
        println(@sprintf "recovery_time = %i,  infection_chance = %.2f,  radius_infection = %.2f" param.recovery_time param.infection_chance param.radius_infection)
        println("flag_multiple_infection = ", flag.flag_multiple_infection ? "true" : "false")
        println("flag_infected_isolation = ", flag.flag_infected_isolation ? "true" : "false")
        println("#================================================================\n")

        out_dir = @sprintf "dat/np=%i_vm=%.3f_vσ=%.3f_vf=%.3f_rii=%i_rt=%i_ic=%.2f_ri=%.2f_fmi=%d_fii=%d/" param.num_particles param.vel_mean param.vel_σ param.vel_flc param.ratio_infection_init param.recovery_time param.infection_chance param.radius_infection flag.flag_multiple_infection flag.flag_infected_isolation
        mkpath(out_dir)
        return out_dir
    end

    """
    Scatter plot of particles
    """
    function plot_particles(itr, param, num_snapshot, out_dir, ptcl)
        x_not_infected = Array{Float64}(undef, num_snapshot.not_infected)
        y_not_infected = Array{Float64}(undef, num_snapshot.not_infected)
        x_infected = Array{Float64}(undef, num_snapshot.infected)
        y_infected = Array{Float64}(undef, num_snapshot.infected)
        x_recovered = Array{Float64}(undef, num_snapshot.recovered)
        y_recovered = Array{Float64}(undef, num_snapshot.recovered)

        # Extract necessary information
        count_not_infeced, count_infected, count_recovered = 1, 1, 1
        for itr_ptcl = 1:param.num_particles
            if ptcl[itr_ptcl].status == "not_infected"
                x_not_infected[count_not_infeced] = ptcl[itr_ptcl].pos_x
                y_not_infected[count_not_infeced] = ptcl[itr_ptcl].pos_y
                count_not_infeced += 1
            elseif ptcl[itr_ptcl].status == "infected"
                x_infected[count_infected] = ptcl[itr_ptcl].pos_x
                y_infected[count_infected] = ptcl[itr_ptcl].pos_y
                count_infected += 1
            elseif ptcl[itr_ptcl].status == "recovered"
                x_recovered[count_recovered] = ptcl[itr_ptcl].pos_x
                y_recovered[count_recovered] = ptcl[itr_ptcl].pos_y
                count_recovered += 1
            end
        end

        itr_str = lpad(itr, 4, "0")
        filename = string(out_dir, "itr_", itr_str, ".png")
        p = scatter(
            x_not_infected, y_not_infected,
            markercolor = :deepskyblue,
            label = "Not infected")
        p! = scatter!(
            x_infected, y_infected,
            markercolor = :orangered,
            label = "Infected")
        p! = scatter!(
            x_recovered, y_recovered,
            aspect_ratio = 1,
            markercolor = :gold,
            label = "Recovered",
            xlims = (0.0, param.x_range),
            ylims = (0.0, param.y_range),
            axis = nothing,
            size=(640, 640),
            title = string("itr = ", itr))
        plot(p)
        # savefig(p, filename)
    end

    """
    Make gif video
    """
    function make_gif(param, anim, out_dir)
        gif(
            anim,
            string(out_dir, "particles.gif"),
            fps=5)
    end

    """
    Plot number of
    - Not infected
    - Infected
    - Recovered
    particle as timeseries
    """
    function plot_num_timeseries(param, num_timeseries, out_dir)
        filename = string(out_dir, "timeseries.png")

        # Array of struct -> array  #####POSSIBLE BETTER SOLUTION?#####
        tseries = Array{Int64}(undef, length(num_timeseries), 3)
        for itr_time = 1:length(num_timeseries)
            tseries[itr_time, 1] = num_timeseries[itr_time].recovered  # Order is altered for visualisation
            tseries[itr_time, 2] = num_timeseries[itr_time].infected
            tseries[itr_time, 3] = num_timeseries[itr_time].not_infected
        end

        p = groupedbar(
            tseries,
            bar_position=:stack,
            bar_width=0.7,
            color = [:gold :orangered :deepskyblue],
            label = ["Recovered" "Infected" "Not_infected"],
            xaxis = ("Time step"),
            yaxis = ("Number of particles"),
            size=(960, 640),
            )
        savefig(p, filename)
    end

end

# ========================================
# Main function
# ========================================

## Declare modules
using ProgressMeter
using Distributions
using Printf
using Plots
gr(
    markerstrokewidth = 0,
    markersize = 10
)
using .ParamVar
using .TimeMarch:
set_initial_condition,
update_particles,
count_status
using .Output:
stdout_condition,
plot_particles,
make_gif,
plot_num_timeseries

# ----------------------------------------
## Declare parameters
# ----------------------------------------
num_particles = 100
max_iteration = 100

x_range = 1.0
y_range = 1.0

vel_mean = 0.01
vel_σ = 0.01
vel_flc = 0.3 * vel_mean

ratio_infection_init = 0.03
recovery_time = 5
infection_chance = 0.3

radius_infection = 0.1

param = ParamVar.Parameters(
    num_particles,max_iteration,
    x_range,y_range,
    vel_mean,vel_σ,vel_flc,
    ratio_infection_init,recovery_time,infection_chance,
    radius_infection)

# ----------------------------------------
## Declare flags
# ----------------------------------------
flag_multiple_infection = true
flag_infected_isolation = false

flag = ParamVar.Flags(
    flag_multiple_infection,
    flag_infected_isolation
)

# ----------------------------------------
## Stdout simulation condition
## & define output directory name
# ----------------------------------------
out_dir = stdout_condition(param, flag)

# ----------------------------------------
## Define array of particle properties
# ----------------------------------------
particles = Array{ParamVar.Particle}(undef, param.num_particles)
for itr_ptcl = 1:param.num_particles
    particles[itr_ptcl] = ParamVar.Particle()
end
# particles .= ParamVar.Particle()  # ERROR: LoadError: MethodError: no method matching length(::Main.ParamVar.Particle)

# ----------------------------------------
## Define array of number of
## not_infected/infected/recovered
## particle in a snapshot & its timeseries
# ----------------------------------------
not_infected, infected, recovered = 0, 0, 0

num_snapshot = ParamVar.NumSnapshot(
    not_infected, infected, recovered)

num_timeseries = Array{ParamVar.NumSnapshot}(undef, 0)  # push! in temporal iteration

# ----------------------------------------
## Set initial condition of particles
# ----------------------------------------
set_initial_condition(param, particles)

# ----------------------------------------
## Time iteration of infection simulation
# ----------------------------------------
progress = Progress(param.max_iteration)
anim = @animate for itr_time = 1:param.max_iteration
    # Update particle properties
    update_particles(param, flag, particles)

    # Count not-infected, infected & recovered number of particles
    num_not_infected, num_infected, num_recovered = count_status(param, particles)

    # Update num_snapshot
    num_snapshot =ParamVar.NumSnapshot(num_not_infected, num_infected, num_recovered)

    # append snapshot values to timeseries
    push!(num_timeseries, num_snapshot)

    # Plot particles for gif video
    plot_particles(itr_time, param, num_snapshot, out_dir, particles)

    # tmp_string = @sprintf "itr_time = %i x[1] = %6.3f y[1] = %6.3f" itr_time particles[1].pos_x particles[1].pos_y
    # println(tmp_string)
    # println("itr_time = ", itr_time, " not infected = ", num_snapshot.not_infected, " infected = ", num_snapshot.infected, " recovered = ", num_snapshot.recovered)
    next!(progress)

    # Finish if there are no infected particles any more
    if num_snapshot.infected == 0
        println("\n No patients at itr_time = ", itr_time, " :Exit")
        break
    end
end

# ----------------------------------------
## Make gif video of simulation
# ----------------------------------------
make_gif(param, anim, out_dir)

# ----------------------------------------
## Make figure of timeseries
# ----------------------------------------
plot_num_timeseries(param, num_timeseries, out_dir)


text:
using Test, YaoBase, LuxurySparse

@testset "batch normalize" begin
    s = rand(3, 4)
    batch_normalize!(s, 1)
    for i = 1:4
        @test sum(s[:, i]) ≈ 1
    end

    s = rand(3, 4)
    ss = batch_normalize(s, 1)
    for i = 1:4
        @test sum(s[:, i]) != 1
        @test sum(ss[:, i]) ≈ 1
    end
end

@testset "hilbertkron" begin
    A,B,C,D = [randn(2,2) for i = 1:4]
    II = IMatrix(2)
    ⊗ = kron
    @test hilbertkron(4, [A, B], [3, 1]) ≈ II ⊗ A ⊗ II ⊗ B
    @test hilbertkron(4, [A ⊗ B, C], [3, 1]) ≈ A ⊗ B ⊗ II ⊗ C
    @test hilbertkron(4, [A ⊗ B], [1]) ≈ II ⊗ II ⊗ A ⊗ B
    @test hilbertkron(4, [A ⊗ B, C ⊗ D], [1, 3]) ≈ C ⊗ D ⊗ A ⊗ B

    U = randn(2,2)
    U2 = randn(4,4)
    m = U2 ⊗ II ⊗ U ⊗ II
    @test m == hilbertkron(5, [U, U2], [2, 4])
end

@testset "batched kron" begin
    A, B, C  = rand(ComplexF64, 4, 4, 3), rand(ComplexF64, 4, 4, 3), rand(ComplexF64, 4, 4, 3)
    D = batched_kron(A, B, C)

    tD = zeros(ComplexF64, 64, 64, 3)
    for k in 1:3
        tD[:, :, k] = kron(A[:, :, k], B[:, :, k], C[:, :, k])
    end

    @test tD ≈ D

    B2 = reshape(transpose(reshape(permutedims(B, (3,1,2)), 3,16)), 4,4,3)
    @test B2 isa Base.ReshapedArray
    @test Array(B2) ≈ B
    D2 = batched_kron(A, B2, C)
    @test tD ≈ D2
end

@testset "random matrices" begin
    mat = rand_unitary(8)
    @test isunitary(mat)
    mat = rand_hermitian(8)
    @test ishermitian(mat)

    @test ishermitian(sprand_hermitian(8, 0.5))
    @test isunitary(sprand_unitary(8, 0.5))
end

@testset "rotmat" begin
    theta = 0.5
    @test rot_mat(ComplexF64, Const.X, theta) ≈ ComplexF64[cos(theta/2) -im * sin(theta/2); -im * sin(theta/2) cos(theta/2)]
    @test rot_mat(ComplexF64, Const.X, theta) |> eltype == ComplexF64
    @test rot_mat(ComplexF32, Const.X, theta) |> eltype == ComplexF32
end


text:
struct DampedSystem{T}
    M::AbstractMatrix{T} # Mass matrix
    K::AbstractMatrix{T} # Stiffness matrix
    D::AbstractMatrix{T} # Internal damping matrix
end

using SparseArrays
function ChangeOfBasis(S::DampedSystem)
    T=Complex{Float64}
    n=size(S.K,1)
    # Solve the GEVD
    Kf=deepcopy(S.K)
    for i=1:n,j=1:n
        Kf[j,i]/=(sqrt(S.M.diag[j]*S.M.diag[i]))
    end
    Ω2,q=eigen(Kf)
    Φ=q./sqrt.(S.M.diag)
    Ω=sqrt.(Ω2)
    # Γ is diagonal in the Φ basis
    # Γ=diag(Φ'*D_int*Φ)
    α=0.02
    Γ=α*Ω
    # Form the matrices Ξ and Q
    Ξ=zeros(T,2n)
    Q=spzeros(T,2n,2n)
    # Solve small 2x2 hyperbolic GEVDs
    # Q is kept as a sparse matrix
    for k=1:n
        Ξ[[k,n+k]],Q[[k,n+k],[k,n+k]]=eigen([0 Ω[k];  Ω[k] Γ[k]], [1.0 0; 0 -1])
        # Transform Q such that Q^T*[I 0;0 -I]*Q=I
        Θ=diag(transpose(Q[[k,n+k],[k,n+k]])*[1.0 0; 0 -1]*Q[[k,n+k],[k,n+k]])
        Q[[k,n+k],k]./=sqrt(Θ[1])
        Q[[k,n+k],n+k]./=sqrt(Θ[2])
    end
    Φ=Φ*Q[n+1:2n,:]
    return Ω,Φ,Ξ,Q
end


text:
"""
    struct TextDiagnosticPrinter <: AbstractDiagnosticConsumer
Hold a pointer to a `clang::TextDiagnosticPrinter` object.
"""
struct TextDiagnosticPrinter <: AbstractDiagnosticConsumer
    ptr::CXDiagnosticConsumer
end

Base.unsafe_convert(::Type{CXDiagnosticConsumer}, x::TextDiagnosticPrinter) = x.ptr
Base.cconvert(::Type{CXDiagnosticConsumer}, x::TextDiagnosticPrinter) = x


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LinearPath

Traverse a spatial object with `N` points in order `1:N`.
"""
struct LinearPath <: AbstractPath end

traverse(object, path::LinearPath) = 1:npoints(object)


text:
# GFF3 Reader
# ===========

mutable struct Reader{S <: TranscodingStream} <: BioGenerics.IO.AbstractReader
    state::BioGenerics.Automa.State{S}
    index::Union{Indexes.Tabix, Nothing}
    save_directives::Bool
    targets::Vector{Symbol}
    found_fasta::Bool
    directives::Vector{Record}
    directive_count::Int
    preceding_directive_count::Int

    function Reader(input::S, index=nothing, save_directives::Bool=false, skip_features::Bool=false, skip_directives::Bool=true, skip_comments::Bool=true) where S <: TranscodingStream

        if isa(index, Indexes.Tabix) && !isa(input.stream, BGZFStreams.BGZFStream)
            throw(ArgumentError("not a BGZF stream"))
        end
        targets = Symbol[]
        if !skip_features
            push!(targets, :feature)
        end
        if !skip_directives
            push!(targets, :directive)
        end
        if !skip_comments
            push!(targets, :comment)
        end
        return new{S}(BioGenerics.Automa.State(input, body_machine.start_state, 1, false), index, save_directives, targets, false, Record[], 0, 0)
    end
end

"""
    GFF3.Reader(input::IO;
                index=nothing,
                save_directives::Bool=false,
                skip_features::Bool=false,
                skip_directives::Bool=true,
                skip_comments::Bool=true)

    GFF3.Reader(input::AbstractString;
                index=:auto,
                save_directives::Bool=false,
                skip_features::Bool=false,
                skip_directives::Bool=true,
                skip_comments::Bool=true)

Create a reader for data in GFF3 format.

The first argument specifies the data source. When it is a filepath that ends with *.bgz*, it is considered to be block compression file format (BGZF) and the function will try to find a tabix index file (<filename>.tbi) and read it if any.
See <http://www.htslib.org/doc/tabix.html> for bgzip and tabix tools.

Arguments
---------
- `input`: data source (`IO` object or filepath)
- `index`: path to a tabix file
- `save_directives`: flag to save directive records (which can be accessed with `GFF3.directives`)
- `skip_features`: flag to skip feature records
- `skip_directives`: flag to skip directive records
- `skip_comments`:  flag to skip comment records
"""
function Reader(input::IO; index=nothing, save_directives::Bool=false, skip_features::Bool=false, skip_directives::Bool=true, skip_comments::Bool=true)

    if isa(index, AbstractString)
        index = Indexes.Tabix(index)
    end

    if isa(input, TranscodingStream)
        return Reader(input, index, save_directives, skip_features, skip_directives, skip_comments)
    end

    stream = TranscodingStreams.NoopStream(input)

    return Reader(stream, index, save_directives, skip_features, skip_directives, skip_comments)

end

function Reader(filepath::AbstractString; index=:auto, save_directives::Bool=false, skip_features::Bool=false, skip_directives::Bool=true, skip_comments::Bool=true)
    if isa(index, Symbol) && index != :auto
        throw(ArgumentError("invalid index argument: ':$(index)'"))
    end
    if endswith(filepath, ".bgz")
        input = BGZFStreams.BGZFStream(filepath)
        if index == :auto
            index = Indexes.findtabix(filepath)
        end
    else
        input = open(filepath)
    end
    return Reader(input, index=index, save_directives=save_directives, skip_features=skip_features, skip_directives=skip_directives, skip_comments=skip_comments)
end

function Base.eltype(::Type{<:Reader})
    return Record
end

function BioGenerics.IO.stream(reader::Reader)
    return reader.state.stream
end

function Base.eof(reader::Reader)
    return reader.state.filled || eof(reader.state.stream)
end

function Base.close(reader::Reader)
    # make trailing directives accessable
    reader.directive_count = reader.preceding_directive_count
    reader.preceding_directive_count = 0
    close(BioGenerics.IO.stream(reader))
end

function Base.read!(reader::Reader, record::Record)
    return readrecord!(reader.state.stream, reader, record)
end

function index!(record::Record)
    stream = TranscodingStreams.NoopStream(IOBuffer(record.data))
    return index!(stream, record)
end

function Base.iterate(reader::Reader, nextone::Record = Record())
    if BioGenerics.IO.tryread!(reader, nextone) === nothing
        return nothing
    end
    return copy(nextone), empty!(nextone)
end

function IntervalCollection(reader::Reader)
    intervals = collect(Interval{Record}, reader)
    return IntervalCollection(intervals, true)
end

function GenomicFeatures.eachoverlap(reader::Reader, interval::Interval)
    if reader.index === nothing
        throw(ArgumentError("index is null"))
    end
    return Indexes.TabixOverlapIterator(reader, interval)
end

"""
Return all directives that preceded the last GFF entry parsed as an array of strings.

Directives at the end of the file can be accessed by calling `close(reader)` and then `directives(reader)`.
"""
function directives(reader::Reader)
    ret = String[]
    for i in lastindex(reader.directives)-reader.directive_count+1:lastindex(reader.directives)
        push!(ret, content(reader.directives[i]))
    end
    return ret
end

"""
Return true if the GFF3 stream is at its end and there is trailing FASTA data.
"""
function hasfasta(reader::Reader)
    if eof(reader)
        return reader.found_fasta
    end
    error("GFF3 file must be read until the end before any FASTA sequences can be accessed")
end

"""
Return a FASTA.Reader initialized to parse trailing FASTA data.

Throws an exception if there is no trailing FASTA, which can be checked using `hasfasta`.
"""
function getfasta(reader::Reader) #TODO: move responsibility to FASTX.jl.
    if !hasfasta(reader)
        error("GFF3 file has no FASTA data")
    end
    return FASTA.Reader(reader.state.stream)
end

function appendfrom!(dst, dpos, src, spos, n)
    if length(dst) < dpos + n - 1
        resize!(dst, dpos + n - 1)
    end
    unsafe_copyto!(dst, dpos, src, spos, n)
    return dst
end

const record_machine, body_machine = (function ()
    cat = Automa.RegExp.cat
    rep = Automa.RegExp.rep
    rep1 = Automa.RegExp.rep1
    alt = Automa.RegExp.alt
    opt = Automa.RegExp.opt

    feature = let
        seqid = re"[a-zA-Z0-9.:^*$@!+_?\-|%]+"
        seqid.actions[:enter] = [:pos]
        seqid.actions[:exit]  = [:feature_seqid]

        source = re"[ -~]+"
        source.actions[:enter] = [:pos]
        source.actions[:exit]  = [:feature_source]

        type_ = re"[ -~]+"
        type_.actions[:enter] = [:pos]
        type_.actions[:exit]  = [:feature_type_]

        start = re"[0-9]+|\."
        start.actions[:enter] = [:pos]
        start.actions[:exit]  = [:feature_start]

        end_ = re"[0-9]+|\."
        end_.actions[:enter] = [:pos]
        end_.actions[:exit]  = [:feature_end_]

        score = re"[ -~]*[0-9][ -~]*|\."
        score.actions[:enter] = [:pos]
        score.actions[:exit]  = [:feature_score]

        strand = re"[+\-?]|\."
        strand.actions[:enter] = [:feature_strand]

        phase = re"[0-2]|\."
        phase.actions[:enter] = [:feature_phase]

        attributes = let
            char = re"[^=;,\t\r\n]"
            key = rep1(char)
            key.actions[:enter] = [:pos]
            key.actions[:exit]  = [:feature_attribute_key]
            val = rep(char)
            attr = cat(key, '=', val, rep(cat(',', val)))

            cat(rep(cat(attr, ';')), opt(attr))
        end

        cat(seqid,  '\t',
            source, '\t',
            type_,  '\t',
            start,  '\t',
            end_,   '\t',
            score,  '\t',
            strand, '\t',
            phase,  '\t',
            attributes)
    end
    feature.actions[:exit] = [:feature]

    directive = re"##[^\r\n]*"
    directive.actions[:exit] = [:directive]

    comment = re"#([^#\r\n][^\r\n]*)?"
    comment.actions[:exit] = [:comment]

    record = alt(feature, directive, comment)
    record.actions[:enter] = [:mark]
    record.actions[:exit]  = [:record]

    blank = re"[ \t]*"

    newline = let
        lf = re"\n"
        lf.actions[:enter] = [:countline]

        cat(opt('\r'), lf)
    end

    body = rep(cat(alt(record, blank), newline))
    body.actions[:exit] = [:body]

    # look-ahead of the beginning of FASTA
    body′ = cat(body, opt('>'))

    return map(Automa.compile, (record, body′))
end)()

const record_actions = Dict(
    :mark => :(@mark),
    :pos => :(pos = @relpos(p)),
    :feature_seqid   => :(record.seqid  = pos:@relpos(p-1)),
    :feature_source  => :(record.source = pos:@relpos(p-1)),
    :feature_type_   => :(record.type_  = pos:@relpos(p-1)),
    :feature_start   => :(record.start  = pos:@relpos(p-1)),
    :feature_end_    => :(record.end_   = pos:@relpos(p-1)),
    :feature_score   => :(record.score  = pos:@relpos(p-1)),
    :feature_strand  => :(record.strand = @relpos(p)),
    :feature_phase   => :(record.phase  = @relpos(p)),
    :feature_attribute_key => :(push!(record.attribute_keys, pos:@relpos(p-1))),
    :feature         => :(record.kind = :feature),
    :directive       => :(record.kind = :directive),
    :comment         => :(record.kind = :comment),
    :record          => quote
        appendfrom!(record.data, 1, data, @markpos, p-@markpos)
        record.filled = 1:(p-@markpos)
    end
)

context = Automa.CodeGenContext(
    generator = :goto,
    checkbounds = false,
    loopunroll = 0
)

Automa.Stream.generate_reader(
    :index!,
    record_machine,
    arguments = (:(record::Record),),
    actions = record_actions,
    context = context,
    returncode = quote
        if cs == 0
            return record
        end
        throw(ArgumentError(string("failed to index ", eltype(record), " ~>", repr(String(data[p:min(p+7,p_end)])))))
    end
) |> eval

Automa.Stream.generate_reader(
    :readrecord!,
    body_machine,
    arguments = (:(reader::Reader), :(record::Record)),
    actions = merge(record_actions,
        Dict(
            :countline => :(linenum += 1),
            :record => quote
                appendfrom!(record.data, 1, data, @markpos, p-@markpos)
                record.filled = 1:(p-@markpos)
                if isfeature(record)
                    reader.directive_count = reader.preceding_directive_count
                    reader.preceding_directive_count = 0
                elseif isdirective(record)
                    reader.preceding_directive_count += 1
                    if reader.save_directives
                        push!(reader.directives, copy(record))
                    end
                    if is_fasta_directive(record)
                        reader.found_fasta = true
                        reader.state.filled = true
                    end
                end

                if record.kind ∈ reader.targets
                    found_record = true
                    @escape
                end
            end,
            :body => quote
                if data[p] == UInt8('>')
                    reader.found_fasta = true
                    reader.state.filled = true
                    # HACK: any better way?
                    cs = 0
                    @goto exit
                end
            end,
        )
    ),
    context = context,
    initcode = quote
        cs = reader.state.state
        linenum = reader.state.linenum
        found_record = false
    end,
    loopcode = quote
        if found_record
            @goto __return__
        end
    end,
    returncode = quote

        reader.state.state = cs
        # reader.state.filled |= cs == 0 # Note: if set to true, remains true.
        reader.state.linenum = linenum

        if found_record
            return record
        end

        if cs == 0 || eof(stream)
            throw(EOFError())
        end

        if cs < 0
            error(eltype(Reader), " file format error on line ", linenum, " ~>", repr(String(data[p:min(p+7,p_end)])))
        end

        if p > p_eof ≥ 0
            error("incomplete $(typeof(reader)) input on line ", linenum)
        end

    end
) |> eval


text:
# Compare simulations with measurements from Schymanski et al. (2017) extracted from this
# Github [repository](https://github.com/schymans/Schymanski_leaf-scale_2016).

using PlantBiophysics
using DataFrames
using CSV
using Plots

# Some constants used in the experiment:
maxiter = 20 # Maximum number of iterations for the algorithm to converge
aₛᵥ = 1 # number of sides used for transpiration (hypostomatous: 1, amphistomatous: 2)
cst = Constants(Cₚ = 1004.834)
A = 20.0 # Not used.

# Defining key functions:
function read_dict(file)
    lines = readlines(file)
    params = Dict{String,Float64}()
    for i in 1:length(lines)
        param = lstrip.(rstrip.(split(lines[i],"=")))
        push!(params, param[1] => parse(Float64,param[2]))
    end
    params
end

function run_simulation!(data,params,aₛᵥ)
    data[!, :Tₗ] .= data[!, :λE] .= data[!, :H] .= data[!, :rbh] .= 0.0
    data[!, :Rn] .= 0.0

    for i in 1:size(data,1)
        meteo = Atmosphere(T = Float64(data.T_a[i]) - params["T0"],
                            Wind = Float64(data.v_w[i]),
                            P = data.P_a[i] / 1000,
                            Rh = data.rh[i])
        leaf = LeafModels(energy = Monteith(aₛᵥ = aₛᵥ, maxiter = maxiter),
                    photosynthesis = ConstantA(A),
                    stomatal_conductance = ConstantGs(0.0, gsw_to_gsc(ms_to_mol(data.g_sw[i],data.T_a[i] - params["T0"],data.P_a[i]/1000))),
                    Rₛ = data.Rn_leaf[i], skyFraction = 2.0, d = data.L_l[i])
        out = energy_balance(leaf,meteo,cst)

        data.Tₗ[i] = out.Tₗ
        data.λE[i] = out.λE
        data.H[i] = out.H
        data.Rn[i] = out.Rn
        data.rbh[i] = 1/out.Gbₕ
    end
    data
end

### Figure 6 a of the article:

# Import the inputs for simulation:
results1_6a = DataFrame(CSV.File("data/schymanski_et_al_2017/results1_6a.csv"))
results1_6a.rh = rh_from_e.(results1_6a.P_wa ./ 1000.0, results1_6a.T_a .+ cst.K₀)
sort!(results1_6a, [:v_w])

# Reading the parameters:
params = read_dict("data/schymanski_et_al_2017/vdict_6a.txt")

# Running the simulation:
run_simulation!(results1_6a,params,aₛᵥ)

scatter(results1_6a.v_w, results1_6a.Elmeas, ylim = (-400,400),
        ylab = "Energy flux from leaf (W m-2)",legend=:inline,
        xlab = "Wind speed (m s-1)", label = "LE meas", color = "blue")
scatter!(results1_6a.v_w, results1_6a.Hlmeas,label = "H meas", color = "red")
scatter!(results1_6a.v_w, results1_6a.Rn_leaf,label = "Rn meas", color = "green")
scatter!(results1_6a.v_w, results1_6a.Hlmeas + results1_6a.Elmeas,label = "H+LE meas", color = "green", shape = :star5)
plot!(results1_6a.v_w, results1_6a.H,label = "H sim", color = "red")
plot!(results1_6a.v_w, results1_6a.λE,label = "LE sim", color = "blue")
plot!(results1_6a.v_w, results1_6a.Rn,label = "Rn sim", color = "green")
savefig("schymanski_et_al_2017_6a.svg")


text:
mutable struct fmpr
  _man::Int
  _exp::Int

  function fmpr()
    z = new()
    ccall((:fmpr_init, :libarb), Void, (Ptr{fmpr}, ), &z)
    finalizer(z, _fmpr_clear_fn)
    return z
  end

  function fmpr(x::Ptr{arf_struct})
    z = new()
    ccall((:fmpr_init, :libarb), Void, (Ptr{fmpr}, ), &z)
    ccall((:arf_get_fmpr, :libarb), Void, (Ptr{fmpr}, Ptr{arf_struct}), &z, x)
    finalizer(z, _fmpr_clear_fn)
    return z
  end

  function fmpr(x::arf_struct)
    z = new()
    ccall((:fmpr_init, :libarb), Void, (Ptr{fmpr}, ), &z)
    ccall((:arf_get_fmpr, :libarb), Void, (Ptr{fmpr}, Ptr{arf_struct}), &z, &x)
    finalizer(z, _fmpr_clear_fn)
    return z
  end
end

function _fmpr_clear_fn(x::fmpr)
  ccall((:fmpr_clear, :libarb), Void, (Ptr{fmpr}, ), &x)
end

function fmpq(x::fmpr)
  z = fmpq()
  ccall((:fmpr_get_fmpq, :libarb), Void, (Ptr{fmpq}, Ptr{fmpr}), &z, &x)
  return z
end

mutable struct cfrac
  coeff::Ptr{fmpz}
  n::Int
  l::Int # real length

  function cfrac(x::Int)
    z = new()
    z.coeff = ccall((:_fmpz_vec_init, :libflint), Ptr{fmpz}, (Int, ), x)
    z.n = x
    finalizer(z, _cfrac_clear_fn)
    return z
  end
end

function _cfrac_clear_fn(x::cfrac)
  ccall((:_fmpz_vec_clear, :libflint), Void, (Ptr{fmpz}, Int), x.coeff, x.n)
end

function show(io::IO, x::cfrac)
  print(io, "[ ")
  for i in 1:x.l
    print(io, unsafe_load(x.coeff, i))
    print(io, " ,")
  end
  print(io, "]")
end

# THIS LEAKS MEMORY
function cfrac(x::fmpq, y::Int)
  r = fmpq()
  z = cfrac(y)
  l = ccall((:fmpq_get_cfrac, :libflint), Int, (Ptr{fmpz}, Ptr{fmpq}, Ptr{fmpq}, Int), z.coeff, &r, &x, y)
  z.l = l
  return z, r
end

function fmpq(x::cfrac)
  z = fmpq()
  ccall((:fmpq_set_cfrac, :libflint), Void, (Ptr{fmpq}, Ptr{fmpz}, Int), &z, x.coeff, x.l)
  return z
end

function fmpq(x::cfrac, y::Int)
  z = fmpq()
  ccall((:fmpq_set_cfrac, :libflint), Void, (Ptr{fmpq}, Ptr{fmpz}, Int), &z, x.coeff, y)
  return z
end


text:
using TransposedVectors
using Base.Test

include("TransposedVector.jl")
include("hcat.jl")
include("mul.jl")
include("ldiv.jl")
include("rdiv.jl")
include("diagonal.jl")
include("triangular.jl")
include("sparse.jl")
include("qr.jl")
include("norm.jl")


text:
"""
Arrays with fixed size and circular indexing.
"""
module CircularArrays

export CircularArray, CircularVector

"""
    CircularArray{T, N, A} <: AbstractArray{T, N}

`N`-dimensional array backed by an `AbstractArray{T, N}` of type `A` with fixed size and circular indexing.

    array[index...] == array[mod1.(index, size)...]
"""
struct CircularArray{T, N, A <: AbstractArray{T,N}} <: AbstractArray{T,N}
    data::A
    CircularArray{T,N}(data::A) where A <: AbstractArray{T,N} where {T,N} = new{T,N,A}(data)
end

"""
    CircularVector{T,A} <: AbstractVector{T}

One-dimensional array backed by an `AbstractArray{T, 1}` of type `A` with fixed size and circular indexing.
Alias for [`CircularArray{T,1,A}`](@ref).

    array[index] == array[mod1(index, length)]
"""
const CircularVector{T} = CircularArray{T, 1}

"""
    CircularArray(data)

Create a `CircularArray` wrapping the array `data`.
"""
CircularArray(data::AbstractArray{T,N}) where {T,N} = CircularArray{T,N}(data)
CircularArray{T}(data::AbstractArray{T,N}) where {T,N} = CircularArray{T,N}(data)

"""
    CircularArray(def, size)

Create a `CircularArray` of size `size` filled with value `def`.
"""
CircularArray(def::T, size) where T = CircularArray(fill(def, size))

Base.IndexStyle(::Type{CircularArray{T,N,A}}) where {T,N,A} = IndexCartesian()
Base.IndexStyle(::Type{<:CircularVector}) = IndexLinear()

@inline Base.getindex(arr::CircularArray, i::Int) = @inbounds getindex(arr.data, mod(i, eachindex(IndexLinear(), arr.data)))
@inline Base.getindex(arr::CircularArray{T,N,A}, I::Vararg{<:Int,N}) where {T,N,A} = @inbounds getindex(arr.data, mod.(I, axes(arr.data))...)

@inline Base.setindex!(arr::CircularArray, v, i::Int) = @inbounds setindex!(arr.data, v, mod(i, eachindex(IndexLinear(), arr.data)))
@inline Base.setindex!(arr::CircularArray{T,N,A}, v, I::Vararg{<:Int,N}) where {T,N,A} = @inbounds setindex!(arr.data, v, mod.(I, axes(arr.data))...)

@inline Base.size(arr::CircularArray) = size(arr.data)
@inline Base.resize!(arr::CircularArray, nl::Integer) = resize!(arr.data, nl)
@inline Base.axes(arr::CircularArray) = axes(arr.data)
@inline Base.parent(arr::CircularArray) = arr.data

@inline Base.iterate(arr::CircularArray, i...) = iterate(parent(arr), i...)

@inline Base.in(x, arr::CircularArray) = in(x, parent(arr))
@inline Base.copy(arr::CircularArray) = CircularArray(copy(parent(arr)))

@inline function Base.checkbounds(arr::CircularArray, I...)
    J = Base.to_indices(arr, I)
    length(J) == 1 || length(J) >= ndims(arr) || throw(BoundsError(arr, I))
    nothing
end

@inline _similar(arr::CircularArray, ::Type{T}, dims) where T = CircularArray(similar(parent(arr), T, dims))
@inline Base.similar(arr::CircularArray, ::Type{T}, dims::Tuple{Base.DimOrInd, Vararg{Base.DimOrInd}}) where T = _similar(arr, T, dims)
# Ambiguity resolution with Base
@inline Base.similar(arr::CircularArray, ::Type{T}, dims::Dims) where T = _similar(arr, T, dims)
@inline Base.similar(arr::CircularArray, ::Type{T}, dims::Tuple{Union{Integer, Base.OneTo}, Vararg{Union{Integer, Base.OneTo}}}) where T = _similar(arr, T, dims)

@inline _similar(::Type{CircularArray{T,N,A}}, dims) where {T,N,A} = CircularArray{T,N}(similar(A, dims))
@inline Base.similar(CA::Type{CircularArray{T,N,A}}, dims::Tuple{Base.DimOrInd, Vararg{Base.DimOrInd}}) where {T,N,A} = _similar(CA, dims)
# Ambiguity resolution with Base
@inline Base.similar(CA::Type{CircularArray{T,N,A}}, dims::Dims) where {T,N,A} = _similar(CA, dims)
@inline Base.similar(CA::Type{CircularArray{T,N,A}}, dims::Tuple{Union{Integer, Base.OneTo}, Vararg{Union{Integer, Base.OneTo}}}) where {T,N,A} = _similar(CA, dims)

@inline Broadcast.BroadcastStyle(::Type{CircularArray{T,N,A}}) where {T,N,A} = Broadcast.ArrayStyle{CircularArray{T,N,A}}()
@inline Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{CircularArray{T,N,A}}}, ::Type{ElType}) where {T,N,A,ElType} = CircularArray(similar(convert(Broadcast.Broadcasted{typeof(Broadcast.BroadcastStyle(A))}, bc), ElType))

@inline Base.dataids(arr::CircularArray) = Base.dataids(parent(arr))

function Base.showarg(io::IO, arr::CircularArray, toplevel)
    print(io, ndims(arr) == 1 ? "CircularVector(" : "CircularArray(")
    Base.showarg(io, parent(arr), false)
    print(io, ')')
    # toplevel && print(io, " with eltype ", eltype(arr))
end

"""
    CircularVector(data)

Create a `CircularVector` wrapping the array `data`.
"""
CircularVector(data::AbstractArray{T, 1}) where T = CircularVector{T}(data)

"""
    CircularVector(def, size)

Create a `CircularVector` of size `size` filled with value `def`.
"""
CircularVector(def::T, size::Int) where T = CircularVector{T}(fill(def, size))

function Base.deleteat!(a::CircularVector, i::Integer)
    deleteat!(a.data, mod(i, eachindex(IndexLinear(), a.data)))
    a
end

function Base.deleteat!(a::CircularVector, inds)
    deleteat!(a.data, sort!(unique(map(i -> mod(i, eachindex(IndexLinear(), a.data)), inds))))
    a
end

end


text:
#=
Cost and final layer delta functions.
These are paired to each type of output layer.
=#

################################################################################
#
# LinearOutputLayer
#
################################################################################

########################################
# Quadratic / [Sum of Squares Diff]       # ℝ
########################################

# J = 1/2 ∑ᵐ ∑ᵏ (h - y)²
# h is a row vector and y is a column vector
function quadratic_cost{F<:AbstractFloat}(h::Array{F}, y::Vector{F})
    #length(y) == size(h, 2) || throw(DimensionMismatch("Inconsistent array sizes"))
    s = F(0)
    @simd for i in eachindex(h)
        @inbounds s += abs2(h[i] - y[i])
    end
    return s * F(0.5)
end

########################################
# Linear Delta                            # ℝ^ 1 x num_obs
########################################

# Final layer delta function for regression
function delta_linear!{F<:AbstractFloat}(Δ::Array{F}, h::Array{F}, y::Vector{F},
                                         num_obs::F)

    # Δ⁽ᴸ⁾ = 1/m (h - y)               (recall: A⁽ᴸ⁾ = h)
    @simd for j in eachindex(y)
        @inbounds Δ[j] = (h[j] - y[j]) / num_obs
    end
end

################################################################################
#
# LogisticOutputLayer
#
################################################################################

########################################
# Cross Entropy                           # ℝ
########################################

# J = ∑ᵐ - y * log(h) - (1-y) log(1-h)
# h is a row vector and y is a column vector
function cross_entropy_cost{F<:AbstractFloat, T<:Integer}(h::Array{F}, y::Vector{T})
    #length(y) == size(h, 2) || throw(DimensionMismatch("Inconsistent array sizes"))
    s = F(0)
    @simd for j in eachindex(h)
        @inbounds if y[j] == T(1)
                      s -= log(h[j])
                  else
                      s -= log1p(-h[j])
                  end
    end
    s
end

########################################
# Logistic Delta                          # ℝ^ 1 x num_obs
########################################

# Final layer delta function for regression
function delta_logistic!{F<:AbstractFloat, T<:Integer}(Δ::Array{F}, h::Array{F},
                                                       y::Vector{T}, num_obs::F)

    # Δ⁽ᴸ⁾ = h - y               (recall: A⁽ᴸ⁾ = h)
    @simd for j in eachindex(y)
        @inbounds Δ[j] = (h[j] - y[j]) / num_obs
    end
end

################################################################################
#
# SoftmaxOutputLayer
#
################################################################################

########################################
# Log Loss                                # ℝ
########################################

# J = ∑ᵐ ∑ᵏ - 1{y=k} log(h)
function log_loss_cost{F<:AbstractFloat, T<:Integer}(H::Array{F}, y::Vector{T})
    length(y) == size(H, 2) || throw(DimensionMismatch("Inconsistent array sizes"))
    s = F(0)
    for j in eachindex(y)
        @inbounds s -= log(H[y[j], j])
    end
    s
end

########################################
# Softmax Delta                          # ℝ^ num_classes x num_obs
########################################

# Final layer delta function for classification
function delta_softmax!{F<:AbstractFloat, T<:Integer}(Δ::Array{F}, H::Array{F},
                                                      y::Vector{T}, num_obs::F)

    # Δ⁽ᴸ⁾ = 1/m (H - Y)               (recall: A⁽ᴸ⁾ = H)
    copy!(Δ, H)
    for j in eachindex(y)
        @inbounds Δ[y[j], j] -= F(1)
    end

    # Divide by m
    broadcast!(/, Δ, Δ, num_obs)
end

################################################################################
#
# MultiLinearOutputLayer
#
################################################################################

########################################
# Quadratic                              # ℝ
########################################

# J = 1/2 ∑ᵐ ∑ᵏ (h - y)²
function quadratic_cost{F<:AbstractFloat}(H::Array{F}, Y::Array{F})
    size(H) == size(Y) || throw(DimensionMismatch("Inconsistent array sizes"))
    s = F(0)
    @simd for i in eachindex(H)
        @inbounds s += abs2(H[i] - Y[i])
    end
    return s * F(0.5)
end

########################################
# Multi Linear Delta                     # ℝ^ num_outputs x num_obs
########################################

# Final layer delta function for regression
function delta_linear!{F<:AbstractFloat}(Δ::Array{F}, H::Array{F}, Y::Array{F},
                                         num_obs::F)

    # Δ⁽ᴸ⁾ = 1/m (h - y)               (recall: A⁽ᴸ⁾ = h)

    @simd for j in eachindex(Y)
        @inbounds Δ[j] = (H[j] - Y[j]) / num_obs
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LLVM_assert_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LLVM_assert")
JLLWrappers.@generate_main_file("LLVM_assert", UUID("1a4ff6f1-5156-5f9a-8aa5-a2840611f610"))
end  # module LLVM_assert_jll


text:
#=
1. Create a regular $N^3$ grid T of random tetrahedra:
	(a) Generate four vertices for each tetrahedron randomly in a half-open fixed-size cube.
	(b) Let these cubes form a regular $N^3$ grid.
2. Create a regular $(N − 1)^3$ grid C of such cubes.
3. Align T and C such that the grid nodes of C are at the centers of the grid cells of T .
4. Measure time for T ∪ C.
=#

using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation;

function t1t2(i,j,k)
	v1 = rand(3,4).+ [2i,2j,2k]
	v2 = rand(3,4).+ [2i,2j,2k]
	ev = [[1,2],[1,3],[1,4],[2,3],[2,4],[3,4]]
	fv = [[2,3,4],[1,3,4],[1,2,4],[1,2,3]]
	t1 = [v1,fv,ev]
	t2 = [v2,fv,ev]
	return t1,t2
end

N = 3  # N ≥ 2
tets = []
for i=0:N-1, j=0:N-1, k=0:N-1
	push!(tets, t1t2(i,j,k)...)
end
V,FV,EV = Lar.struct2lar(Lar.Struct(tets))
GL.VIEW([ GL.GLGrid(V,FV, GL.COLORS[1],1), GL.GLFrame2 ]);

W,(_,EW,FW,_) = Lar.cuboidGrid([N-1,N-1,N-1],true);
U = (W .* 2) .+ 0.5
GL.VIEW([ GL.GLGrid(V,FV,GL.COLORS[1],0.5), GL.GLGrid(U,FW,GL.COLORS[1],0.5), GL.GLFrame2 ]);

V,FV,EV = Lar.struct2lar(Lar.Struct([Lar.Struct(tets), Lar.Struct([(U,FW,EW)]) ]))
GL.VIEW([ GL.GLGrid(V,EV, GL.COLORS[1],0.5), GL.GLFrame2 ]);

	cop_EV = convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells));
	cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
	W = convert(Lar.Points, V');
	V, copEV, copFE, copCF = Lar.space_arrangement( W, cop_EV, cop_FE);
	W = convert(Lar.Points, V');
	V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF);

GL.VIEW( GL.GLExplode(V,CVs[1:end],5,5,5,99,0.5) );
GL.VIEW( GL.GLExplode(V,CVs[4:4],5,5,5,99,0.5) );


text:
using KernelFunctions
using MLDataUtils
using Zygote
using Flux
using Distributions, LinearAlgebra
using Plots

Flux.@functor SqExponentialKernel
Flux.@functor KernelSum
Flux.@functor Matern32Kernel
Flux.@functor FunctionTransform

neuralnet = Chain(Dense(1, 3), Dense(3, 2))
k = SqExponentialKernel(FunctionTransform(neuralnet))
xmin = -3;
xmax = 3;
x = range(xmin, xmax; length=100)
x_test = rand(Uniform(xmin, xmax), 200)
x, y = noisy_function(sinc, x; noise=0.1)
X = reshape(x, :, 1)
λ = [0.1]
function f(x, k, λ)
    return kernelmatrix(k, X, x; obsdim=1) *
           inv(kernelmatrix(k, X; obsdim=1) + exp(λ[1]) * I) *
           y
end
f(X, k, 1.0)
loss(k, λ) = ŷ -> sum(y - ŷ) / length(y) + exp(λ[1]) * norm(ŷ)(f(X, k, λ))
loss(k, λ)
ps = Flux.params(k)
# push!(ps,λ)
opt = Flux.Momentum(1.0)
##
for i in 1:10
    grads = Zygote.gradient(() -> loss(k, λ), ps)
    Flux.Optimise.update!(opt, ps, grads)
    p = Plots.scatter(x, y; lab="data", title="Loss = $(loss(k,λ))")
    Plots.plot!(x, f(X, k, λ); lab="Prediction", lw=3.0)
    display(p)
end


text: rm("crud.txt")

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming


text:
using DelegationPattern
using Test

@testset "DelegationPattern.jl" begin
    include("funcs.jl")

    # Note: these functions just output test results to console
    @test test_account() == nothing
    @test test_savings_account() == nothing
    @test test_savings_account2() == nothing
end


text:
using Boltzmann
using MNIST

X, y = traindata()
X = X[:, 1:1000]                     # take only 1000 observations for speed
X = X / (maximum(X) - (minimum(X)))  # normalize to [0..1]

layers = [("vis", GRBM(784, 256)),
          ("hid1", BernoulliRBM(256, 100)),
          ("hid2", BernoulliRBM(100, 100))]
dbn = DBN(layers)
fit(dbn, X)
transform(dbn, X)

dae = unroll(dbn)
transform(dae, X)

save_params("test.hdf5", dbn)
save_params("test2.hdf5", dae)
rm("test.hdf5")
rm("test2.hdf5")


text:
@testset "ppi" begin
    g, train_X, train_y, train_ids = traindata(PPI())
    @test typeof(g) == SimpleDiGraph{Int32}
    @test nv(g) == 44906
    @test ne(g) == 1271267
    @test typeof(train_X) == Array{Float32,2}
    @test size(train_X) == (44906, 50)
    @test typeof(train_y) == SparseMatrixCSC{Int32,Int64}
    @test size(train_y) == (44906, 121)

    g, valid_X, valid_y, valid_ids = validdata(PPI())
    @test typeof(g) == SimpleDiGraph{Int32}
    @test nv(g) == 6514
    @test ne(g) == 205395
    @test typeof(valid_X) == Array{Float32,2}
    @test size(valid_X) == (6514, 50)
    @test typeof(valid_y) == SparseMatrixCSC{Int32,Int64}
    @test size(valid_y) == (6514, 121)

    g, test_X, test_y, test_ids = testdata(PPI())
    @test typeof(g) == SimpleDiGraph{Int32}
    @test nv(g) == 5524
    @test ne(g) == 167461
    @test typeof(test_X) == Array{Float32,2}
    @test size(test_X) == (5524, 50)
    @test typeof(test_y) == SparseMatrixCSC{Int32,Int64}
    @test size(test_y) == (5524, 121)
end


text: using Test, ComputerVisionWithJuliaBook

max_stars_repo_name: bemcho/ComputerVisionWithJuliaBook

max_issues_repo_name: bemcho/ComputerVisionWithJuliaBook

max_forks_repo_name: bemcho/ComputerVisionWithJuliaBook